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Abstract 

Current and future astronomical survey facilities provide a remarkably rich opportunity for transient astronomy, 
combining unprecedented fields of view with high sensitivity and the ability to access previously unexplored wavelength 
regimes. This is particularly true of LOFAR, a recently-commissioned, low-frequency radio interferometer, based in the 
Netherlands and with stations across Europe. The identification of and response to transients is one of LOFAR’s key 
science goals. However, the large data volumes which LOFAR produces, combined with the scientific requirement for 
rapid response, make automation essential. To support this, we have developed the LOFAR Transients Pipeline, or 
TraP. The TraP ingests multi-frequency image data from LOFAR or other instruments and searches it for transients and 
variables, providing automatic alerts of significant detections and populating a lightcurve database for further analysis 
by astronomers. Here, we discuss the scientific goals of the TraP and how it has been designed to meet them. We 
describe its implementation, including both the algorithms adopted to maximize performance as well as the development 
methodology used to ensure it is robust and reliable, particularly in the presence of artefacts typical of radio astronomy 
imaging. Finally, we report on a series of tests of the pipeline carried out using simulated LOFAR observations with a 
known population of transients. 

Keywords: astronomical transients, time domain astrophysics, techniques: image processing, methods: data analysis, 
astronomical databases 


1. Introduction 

1.1. Slow transients and variable sources 

While most objects in the Universe are steady on human 
timescales, there are classes of sources displaying variability 
on timescales of years to days, seconds or fractions of a 
second. This variability can occur regularly, at irregular 
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intervals, or as a singular event from a given object. Search¬ 
ing for these transients and variables requires observatories 
with a large field of view, a capability which was up to 
now reserved only for some optical telescopes and X- and 
q-ray satellites. The radio regime is now also entering this 
area of time-domain astronomy, with several new facilities 
being built that have large fields of view (several square de¬ 
grees or larger) and transients as one of their key scientific 


objectives (e.g., |Taylor et ah 2012 Murphy et al., 2013 
~]van Haarlem et al. 1 12013[ |Bell et al. 


|Tingay et al.[ |2013 
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2014|. A few of these observatories are probing the low 


radio frequency regime, from tens to hundreds of MHz, a 
range that has been largely unexplored so far. 

Transients in the low-frequency radio sky can be divided 
into roughly two classes, characterized by their emission 
processes and the observing techniques used to study them: 
coherent emitters, which display very fast variability, and 
are found mostly in beamformed time series data, and 
incoherent emitters, which display slow variability and are 
usually detected by comparing multiple images of the same 


1.2. Detecting transients and variables 

The transient sky has long been studied across the elec¬ 
tromagnetic spectrum, but the scale of transient searches 
has increased markedly recently, in particular in the optical 
and radio regimes. 

Searching for transients with large field-of-view X- and 
7 -ray instruments has been common for a long time, and a 
variety of techniques have been used for all-sky monitors 


on board the Rossi X-ray Timing Explorer (Levine et al. 


field (Fender and Bell 2011). Here, we take the dividing 
line between slow and fast as ~ 1 second, which is the 
fastest time scale at which radio images are typically made 


1996 

, Compton Gamma Ray Observatory ( 

Fishman 

1992 


Swift (Gehrels et al. 2004), and Fermi Gamma-ray Space 


Telescope (Atwood et al. 2009 Meegan et al. 2009). The 


(but see e.g. Law et al. 2011). The most well known 
examples of coherent emitters are pulsars and masers, but 
coherent emission processes are also predicted, albeit not 
yet discovered, for other sources like gamma-ray bursts 


most common way to find rapid transients at these energies 
is by monitoring a large fraction of the sky, and triggering 
on a sudden increase in the total X- or 7 -ray flux. Alter¬ 
native techniques are required for transients that evolve 
more slowly: for instance, the Earth occultation method 


magnetar flares (Lyubarsky 2014). This paper, however, 


(Usov and Katz 2000 Sagiv and Waxman 2002) and described by Harmon et al. (2002). 


focuses on searching for incoherent transients and variable 
sources in the image plane, on timescales from seconds to 
years. 

The main incoherent emission process at low radio 
frequencies is synchrotron radiation, which arises when 
relativistic electrons are accelerated in strong magnetic 
fields. It is produced where a large amount of energy is in¬ 
jected into the ambient medium in jet sources and explosive 
events, such as X-ray binaries, active galactic nuclei, tidal 
disruption events, gamma-ray bursts, supernovae, mag- 


netars, and flare stars (e.g., 

Dent 1965 

Gregory et al. 

1972 

Frail et al. 

1997, 1999 

Levan et al. 

2011 

). While 


rnany of these sources show short bursts of emission at X- 
or 7 -ray energies, their variability timescale at low radio 
frequencies is much longer, because the radiative lifetimes 
of the particles to synchrotron emission are very long, and 


due to synchrotron self-absorption effects (van der Laan 


1966). Although the latter decreases the sources’ bright¬ 


ness, making their detection more challenging, it has a 
high pay-off scientifically since determining the evolution 
of the spectrum at low radio frequencies provides important 
information on the energetics involved in these events, the 
acceleration of electrons up to relativistic velocities, the gen¬ 
eration of magnetic fields, the production and collimation 
of jets, and the feedback of these jets on their surround¬ 
ings. Furthermore, sources with small angular scales on 
the sky, like active galactic nuclei, show variability which 
is not intrinsic but caused by scattering in the interstellar 


medium (R.ickett 1990). Therefore these sources are not 


only interesting for studying their physical properties, but 
can also be used to probe the medium in between them and 
us. In this context we note that some coherent events that 
are intrinsically very short can be scattered and dispersed 
in the interstellar medium, smearing out their signal to 
timescales that are probed by image transient searches (see 
e.g. 


Broderick et al. 


m prep. 


Transient searches in the image domain over similarly 
large fields-of-view are now planned—and, indeed, already 
being carried out—at optical and radio frequencies. Here, 
efficiently searching the extremely large data volumes pro¬ 
duced is challenging. Optical telescopes optimized to search 
for transients include the Catalina Real-Time Transient 


Survey (I 

Irake et al. 

2009 

), Palomar Transient Factory 

Rau et al 

,2009 

1 , Pan-STARRS ( 

Denneau et al. 

2013 

I, and 


the Las Cumbres Observatory Global Telescope Network 


(Brown et al. 2013). Several radio telescopes have dedi¬ 


cated transient programs as well, notably the Jansky Very 


Large Array (JVLA), AMI (Staley et al. 20131, MeerKAT 


(Karoo Array Telescope; Booth and Jonas 2012), ATA 
(Allen Telescope Array; | Welch et al. 2009) ASKAP (the 
Australian Square Kilometre Array Pathfinder; |Murphy 


et aL||2013), the MWA (Murchison Wideficld Array; Tingay 
eFaLj]20l3|[Bell et al.[|20Mj |, the LWA (Long Wavelength 


Array; Taylor et al. 2012) and LOFAR (the Low Frequency 
Array; van Haarlem et al. 2013). In the longer term, the 


Large Synoptic Survey Telescope (LSST; Ivezic et al. 2014) 


in the optical and Square Kilometre Array (SKA; Dewdney 


et al. 2010 ) will produce a dramatic increase in the number 


of transients which can be detected. 

Broadly, there are two possible techniques which are 


adopted by these searches: difference imaging (e.g. Alard 
and Lupton 1998; Law et ah] 2009) or on a comparison of 


a list of sources measured in a given image against a deep 


reference catalogue (e.g. Drake et al. 2009). 

Difference imaging has been demonstrated to be effec¬ 
tive when applied to optical data, particularly in crowded 
fields which would suffer from source-confusion in a catalogue- 
based survey. However, the efficacy of difference imaging 
in the optical is partly due to to the sources of noise being 
relatively well characterised, with pixel-noise largely inde¬ 
pendently distributed and occurring on a different spatial 
scale to real sources (assuming a well-sampled point-spread 
function), and the fact that optical survey point-spread 
functions usually vary in a smooth fashion amenable to 
model-fitting. 
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In contrast, noise in radio-synthesis images is inher¬ 
ently correlated on similar scales to the sources of interest. 
Furthermore, effects such as radio frequency interference 
(RFI) and interaction between faint beam-sidelobes and 
bright out-of-field sources may cause artefacts which are 
harder to characterise and correct for than those found in 
optical data. As a result, higher signal-to-noise thresholds 
are typically applied to ensure that most spurious detec¬ 
tions are rejected (although this process remains fallible; 


Frail et al. 2012). This degrades the sensitivity advantage 


of the difference imaging technique, and so a cataloguing 
survey provides equivalent results with the added benefit 
of recording lightcurves for individual sources. 

Many recent developments, including the precursors 
of this work, focus on the latter approach: compiling 
lightcurves, storing them in a database, and then search¬ 
ing for transients with a variety of statistical techniques 


(Spreeuw, 2010 Bannister et al.| 2011; Bower et al.J 

2011 

Groff et al. 2011 Swinbank 2011; Thyagarajan et al. 

2011 

Banyer et al. 2012 

Hancock et al. 2012 Williams et al. 

2013 

Mooley et al. 

2013 

Bell et al. 2014 

). The same 


strategy has been adopted in this work, which describes 
the pipeline and methods developed for searching for tran¬ 
sients with LOFAR. The system described here also has a 
broader applicability to other instruments and is developed 
with an eye to the long-term requirements of the SKA. 


2. Transients with LOFAR 

2.1. The International LOFAR Telescope 

LOFAR is a recently-commissioned radio interferometer 
based in the Netherlands and with stations across EuropcQ 
LOFAR operates in the low-frequency radio regime, observ¬ 
ing in the frequency ranges 10-80 (low band) and 110-240 
(high band) MHz, corresponding to wavelengths between 30 
and 1.2 metres. The system pioneers the concept of a soft¬ 
ware telescope , as signals received by simple crossed-dipole 
antennas, which are sensitive to the whole sky, are digitized 
as soon as possible and signal processing is done in software. 
In the low band the voltages from each dipole are directly 
digitized; in the high band, an analogue combination of 
data from 16 antenna elements (a “tile”) is formed prior 
to digitization. 

In a typical imaging observation, each station acts as a 
phased array: the digitized data from each of the antennas 
in the station is coherently combined (“beamformed”) to 
point the station in a particular direction. The field of view 
of beams formed in this way depends on frequency and 
station configuration, with the full width at half maximum 
ranging from 19.55° for a core Dutch LOFAR station at 
30 MHz to 1.29° for an international LOFAR station at 
240 MH^j This beamformed data is then transported to 


the central processing (CEP) facility at the University of 
Groningen where it is correlated by a GPU-based system. 

The software-based processing model provides for a 
great deal of flexibility. After digitization, a polyphase filter 
splits the data into 0.2 MHz wide subbands, the number 
of subbands depending on the quantization of the data: 
it is possible to trade off dynamic range for increased 
subband number and hence bandwidth: in “16 bit mode”, 
244 subbands are available; in “8 bit mode”, 488. These 
subbands can be spread across frequency space, to give 
a large observing bandwidth (48.8MHz in 16 bit mode). 
Alternatively, the beamformer can be configured to form 
multiple beams with different selections of subbands in 
different directions. In this latter mode, by trading off 
against observing bandwidth, an extremely wide field of 
view may be observed. 

When operating in imaging mode, the correlator pro¬ 
vides a dump time of 1 s, and it is this which provides 
a lower-limit to the timescales which can be searched for 
image plane transients. An alternative is to directly search 
high time resolution beamformed data for fast transients, 
as described by |Stappers et al.| (|201ip and |Coenen et al.J 
(2014). It is also possible to configure the telescope such 
that beamformed and image data is recorded simultane¬ 
ously, providing the greatest possible range of cadences in 
a transient search; ultimately, continual monitoring in this 
mode is an important goal. 


2.2. The Transients Key Science Project and the Radio 
Sky Monitor 

LOFAR’s development and commissioning have been 
driven by six science areas: the epoch of reionization, deep 
extragalactic surveys, cosmic magnetism, solar physics and 
space weather, cosmic rays, and transient and variable 
sources. The last of these is the remit of the Tran sients 
Key Science Projeclj^] (TKSP; Fender et al. 2006). The 
TKSP’s interests include transient and variable sources 
on all timescales, from sub-second changes in beamformed 
data (Stappers et al. 2011) to multi-year variability mon¬ 


itored through long-term imaging programmes; see van 


Haarlem et al. (2013) and references therein for a complete 


discussion of the TKSP science case. It is upon detec¬ 
tion and monitoring in the image plane which this work 
concentrates. 

A key programme for the TKSP is the “Radio Sky 
Monitor”, or RSM ( |Fender et al.| 


m prep. 


In this mode, 

multiple beams from LOFAR are used to tile out a large 
area on the sky (Fig. [l]). This field of view is then imaged 
on a logarithmic range of timescales, from 1 to 10000 s, and 
at a range of frequencies, and that image stream is searched 
for transient and variable sources. The survey strategy is 
flexible, but most plausible strategies will focus on the 


1 http://www.astron.nl/-heald/lofarStatusMap.html 
•'The various types of LOTAH stations together with their key pa¬ 
rameters are listed at http://www.astron.nl/radio-observatory/ 


astronomers/lofar-imaging-capabilities-sensitivity/ 
lof ar-imaging-capabilities/lof a see |van Haarlem et all] 

for detailed information on LOFAR’s configuration, 
http: //www. transientskp. org/ 
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Figure 1: The Radio Sky Monitor concept. Multiple LO- 
FAR station beams tile out a large field of view, while other 
beams are available for target observations. 


galactic plane and the zenith, while taking advantage of 
the large field of view to ensure that a large fraction of 
the sky is regularly monitored. While this procedure is 
ongoing, individual beams can be diverted from the ongoing 
survey to monitor specific targets of interest or respond to 
transient alerts in real time (although the latter is currently 
not implemented). 

There are two key data products which result from this 
RSM: an archive of the lightcurves observed for all point 
sources in the LOFAR sky, and low-latency alerts of tran¬ 
sient events. It is the TKSP’s policy that, in general, these 
products will both be made available to the community at 
large. 

While the RSM is running, a large volume of correlated 
(visibility) and image data will be generated. It is regarded 
as impractical to archive all of this data. Instead, an 
averaged version of the visibility data may be stored with 
reduced time and/or frequency resolution, and thumbnail 
images of significant detections recorded. 

Ultimately, LOFAR is designed to provide low-latency 
“streaming” image data. When this is available, the tran¬ 
sient search may be run, and alerts produced, in real time. 
At time of writing, however, this capability is still in devel¬ 
opment. Instead, visibility data is stored to disk for later 
“offline” imaging. This non-real-time version of the system 
has been deployed for LOFAR’s initial operational phase. 
In this mode, visibility data collected by LOFAR undergoes 
some initial post-pro cessing at CEP before being delivered 
to the TKSP. Project members then image this data on 
local hardware, before running the images through a ver¬ 
sion of the transients detection system which is optimized 


for offline use. In this way, TKSP members are able to 
develop our understanding of LOFAR’s imaging capabili¬ 
ties and to test and commission the transients detection 
and monitoring pipeline (or “TraP”) in advance of its full 
deployment as part of a real-time LOFAR system. 

It is on this TraP system which this manuscript focuses. 
In ([3] we provide an overview of its inputs, outputs and 
overall design. In <|4] we describe in detail the algorithms 
employed by the key pipeline components, and in fj5] de¬ 
scribe the data products the pipeline delivers. Section [6] 
describes the pipeline implementation on a technical level. 
Section [7] discusses the development approaches taken. In 
Sj8]we describe testing carried out on the TraP with simu¬ 
lated datasets. Finally, in j 10 we describe enhancements 
which are planned for future releases. 

The TraP was developed with the aim of finding and 
monitoring transients in RSM-like data. However, it is 
worth emphasizing that it should be ultimately applicable 
to a much wider range of instrumentation. For example, it 
is planned to use the TraP to scan as much LOFAR imag¬ 
ing data as possible, in a so-called “piggyback” mode. An 
early version of the TraP has already been used in a study 
of archival VLA data (Bell et al. 2011), while a variant 
will also be deployed as the transient detection system for 
AARTFAAC (the Amsterdam-ASTRON Radio Transients 
Facility and Analysis Centre, an all-visible-sky monitor 
operating commensally with LOFAR; Prasad and Wijiij 


holds 2012). Other developments target the Arcminute 


Microkelvin Imager Large Array (AMI-LA; Staley et al. 


2013 Anderson et al. 20141, a 15 GHz aperture synthesis 


aperture synthesis radio telescope near Cambridge in the 
UK, and KAT-7/MeerKAT, SKA-precursor telescopes in 
the Karoo Desert, South Africa. Further variants targeting 
optical data are also under consideration. 

The TraP is available as open source software; for more 
details, refer to (IT and the code repositorjj/] 


3. Pipeline overview 

The design goal of the TraP is to automatically and 
rapidly identify transient and variable sources within a 
time-series of image data. These sources may be identified 
in two ways: 

• New detections are sources which appear at a location 
where, in previous epochs, no source was seen; 

• Variables are sources which have been observed for 
multiple epochs and show significant variability in 
their lightcurves. 

Such sources are identified automatically by the TraP, 
based on no prior knowledge. It is also possible for the user 
to specify the location of known sources for monitoring. 
Variability metrics are retained for all sources, so that 


J https://github.com/transientskp/tkp/ 
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decisions on what constitutes an ‘interesting source’ may 
be made after data-processing (([ 5 ]). 

Since the TraP is ultimately designed to perform near 
real-time processing of an image stream, we assume that 
after an image has been processed it is no longer available 


for further analysis (modulo the system described in f 6.5). 


Therefore, the TraP currently provides no capability to 
look back at previously processed images in the light of 
new data: it does not, for example, attempt to go back and 
check earlier images of the same position after observing a 
new transient. Although retaining the full image stream is 
unlikely to be practical for projects which generate substan¬ 
tial data volumes, future versions of the TraP may include 
the capability to generate and store an average of the input 
data, using this to increase the depth of the survey and 
improve source characterization. 

3.1. Inputs 

The fundamental input to the trap is a time-series 
of three-dimensional (two spatial, one frequency) image 
“cubes”. These are generally assumed to be produced by the 
LOFAR imaging pipeline ( Heald et al.| 2010 2011| 2014 


van Haarlem et al., 2013), however, as described in (4.1 


the code is designed to be agnostic as to the format and 
origin of the data being ingested. 

In addition, the TraP may optionally be given a user- 
defined list of monitoring positions. Measurements are 
made and stored for each such position in each plane of 
every image cube ingested, regardless of whether the auto¬ 
matic source-finding routines regard it as significant. 


3.2. Products 

The TraP is designed to produce two key data products: 

• Near real-time alerts to the community and/or tar¬ 
geted at specific partners describing ongoing transient 
celestial events; 

• An archival database of lightcurves for all astronomi¬ 
cal point sources detected during pipeline processing 
together with information about their variability. 


The pipeline system is flexible enough to provide alerts 
in a variety of formats, and it is therefore able to interoper¬ 
ate with whatever mechanisms other facilities have in place 
for receiving notifications. For example, one can imagine 
e-mail or SMS being convenient. However, development 


has focused on the VOEvent system (Seaman et al. 20111 


and its association distribution networks (Williams et al. 


2012). These provide a flexible and convenient method for 


widespread alert dissemination, which is described in detail 
in ( ]6.6| 

In addition to these fundamental data products, the 
TraP may optionally store a copy of all the image pixel 
data processed for future reference. This is not required 
for the analysis performed by the TraP, but we have found 
it convenient to maintain an archive of some or all of 
the images processed for display purposes (e.g. using the 
interface described in (5.2.2). 



Figure 2: Transients Pipeline (“TraP”) overview, showing 
the flow of data through the system as an image “cube” 
is processed, ultimately producing an archival lightcurve 
database and (if appropriate) transient alert messages. The 
individual pipeline stages are described in and El their 
implementation, in terms of a mixture of Python code and 
database routines, is discussed in (Jbj The dotted parts of 
the diagram represent optional functionality: they are not 
required for the core TraP functionality. 


3.3. Methods 

To map from the inputs to the products described above, 
the following procedure is adopted. Each of these stages is 
described in more detail in (|4j their relationship is shown 
graphically in Fig. [2] 


1 . 


2 . 


3. 


The image cube is stored in main memory as a se¬ 
ries of two-dimensional images, each representing a 
particular frequency. The in-memory representation 
of an image used by the TraP is independent of the 
on-disk data storage format; see (]4.1|for details. 


Each image undergoes a quality-control procedure, 
designed to identify and exclude from further pro¬ 
cessing data of unacceptably poor quality. Note that 
even when data is excluded by quality control it is 
not completely disregarded, but rather its existence 
and the reason for its rejection are recorded. For 
details of the quality control checks and the way in 


which they are applied see (4.2 


A source-finding and measurement process is ap¬ 
plied independently to each plane of the image cube. 
Sources which meet the relevant significance criteria 
are parameterized by elliptical Gaussians. For more 


details on the source finding procedure see (4.3 
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4. 


5. 


6 . 


An “association” procedure is carried out, in which 
the new measurements are either identified as updates 
to the lightcurves of previously known sources, or as 
new, previously-undetected, sources. Details on the 
algorithms used for source association may be found 
in §4.4| 

A fist of sources which are expected to appear within 
the image but were not detected by the source find¬ 
ing procedure above is now constructed following the 
procedure described in 14.5 The same source mea¬ 


surement code is now used to fit and record source 
parameters at each of these fixed positions, and the 
relevant lightcurves updated accordingly. 

The same source measurement code is used to fit and 
record source parameters at each of the user-specified 
monitoring positions, and the relevant lightcurves 
updated accordingly. This procedure is described in 


after the insertion of any particular image cube. 

Finally, it should be noted that although it is possi¬ 
ble to create LOFAR images with full polarization, and 
notwithstanding the ultimate TraP design goals, the cur¬ 
rent version of the TraP searches for transient and variable 
sources only within total intensity (Stokes I) images, and 
other polarization information is not used. In time, though, 
polarization information will be essential for properly char¬ 
acterizing the sources being identified: see j ]10| for more 
information on future plans. 

4. Key pipeline stages 

In this section we describe the logical structure of the 
TraP, focusing on the core stages of the pipeline and the 
algorithms that they employ. Section [6] describes how 
this logical structure is implemented in terms of deployed 
software and hardware resources. 


7. For each lightcurve in the database a series of aggre¬ 
gate properties describing the astronomical source 
are calculated. These include a weighted mean posi¬ 
tion, flux density and a series of variability indices 
which quantify the magnitude and significance of the 
variability of the lightcurve. This is described in §4.7| 

At this point, the core loop of the pipeline is complete: 
the next (by time) image cube in the sequence may be 
ingested and the process repeats. At the same time, the 
results are immediately made available via a database, as 
described in ^5} 

Further analysis may be performed by querying the 
database for scientifically relevant phenomena and reacting 
appropriately. For example, one could search for all bright 
sources which do not have a previously-detected counter¬ 
part, and thereby identify new transients. Alternatively, 
a query could search for lightcurves which have reached 
a particular threshold in the variability indices, or which 
meet some other user-defined criteria. 

It is important to emphasize that these queries can be 
performed at any time. For example, the user could wait 
until the complete pipeline run has been completed and 
archived before searching the database; equally, however, a 
real-time analysis system can query the database continu¬ 
ously as new results are added, and thereby identify new 
transients immediately. 

As new measurements are appended to the database, 
continuously-updated measures such as the variability in¬ 
dices for a given lightcurve or the weighted mean position 
of the associated astronomical source will change with time. 
It is possible, therefore, that a particular source which was 
identified as variable by the real-time analysis system at 
some particular timestep will, in the fullness of time, be 
shown to not, in fact, vary significantly. In order to ensure 
reproducibility, the database records all the intermediate 
values as well as the final, archival result. That is, the user 
may query the archival database not just for the eventual 
state of a particular source, but for its state as recorded 


4-1- Data accessors 

While the TraP has been developed with LOFAR in 
mind, many of the core issues we are addressing are widely 
applicable to much of the emerging field of transient astron¬ 
omy. As such, we aim to make it easy to adapt the TraP 
to ingest images from data sources other than LOFAR. 
The pipeline is designed to be data-source agnostic: the 
origin of the data is abstracted away from the scientific 
logic. This has proven to be useful as internal LOFAR 
data-storage standards have evolved. 

Data source abstraction is achieved by defining a uni¬ 
form interface which all routines in the pipeline use to access 
and manipulate image data. Derived classes represent data 
from specific sources, providing different routines for load¬ 
ing the data, translating telescope-specific metadata, and 
so on. Adding support for data from a new telescope is 
generally straightforward: for most instruments, just a few 
simple extensions to the predefined routines for working 
with FITS0 or CASA0/ casacor^jare required. 

This system has enabled the TraP to be used in con¬ 
junction with data not only from LOFAR but also from 


the VLA and AMI-LA, as described in (2.2 


4-2. Quality control 

During the imaging procedure and the future real-time 
imaging pipeline, extremely large numbers of images will be 
produced for processing by the TraP. Some of these images 
will be of insufficient standard for transient searches: for 
instance, high RFI or poor calibration solutions can lead 
to increased noise in the image or artefacts that may be 
mistaken as transients. 

The quality control procedure identifies and rejects 
those images which do not pass a range of tests. The sys¬ 
tem is modular: new tests can easily be added as required. 


5 http://fits.gsfc.nasa.gov/ 
e http://casa.nrao.edu/ 

'http://casacore.googlecode.com/ 
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Further, tests may be made specific to certain instrumenta¬ 
tion by building upon the data accessor framework (j |4.1[ ). 

The standard checks supplied with the released version 
of the TraP are all LOFAR specific. They are: 

• Test that the measured noise in the image does not 
significantly exceed the theoretically expected value 

• Test for appropriate sampling and shape of the restor¬ 
ing beam parameters (j |4.2.2 1 ; 

• Test for proximity of the image pointing direction to 
bright radio sources (j |4.2.3 ). 

An image which fails one or more of these tests is not 
further processed. Details of the failure are logged for 
future reference. 

These tests are designed to provide a quick and simple 
mitigation of common failures observed during develop¬ 
ment and commissioning. As TraP moves into production 
deployments, it will be possible to supplement them with 
a range of more elaborate tests as and when required. 

4-2.1. Check for noisy images 

A clear signature of a poor quality image is when the 
measured noise level significantly differs from the theoret¬ 
ically expected value: measured values which are either 
too low or too high are indicative of problems with the 
observation or its calibration. 

The theoretical noise in LOFAR images can be calcu¬ 
lated using parameters extracted from the image metadata, 
such as the array configuration and integration time used 


we conduct 

1. Select a square region congruent with the image cen¬ 
tre and comprising 25 % of the total pixels in the 
image; 

2. Iteratively reject pixel values more than n standard 
deviations from the median, where n is some user- 
defined parameter (typically four), until no further 
pixels are being rejected. 

3. Using the remaining pixels, calculate the mean pixel 
value and the RMS scatter around this mean value. 

We then calculate a simple ratio between the measured 
RMS noise and the theoretical noise. The image is rejected 
when this ratio falls outside a user-specified range. 


(Nijboer et al. 2009 van Haarlem et al. 2013 


To measure the observed RMS in an image, 
the following steps: 


4-2.2. Check restoring beam parameters 

The properties of the restoring beam (Hogbom 19741 
used to create the images used within the TraP also play a 
significant role in assessing the image quality. The image 
should be created such that the beam is appropriately 
sampled, with around three pixels across its minor axis. 
Incorrect sampling can cause increased noise, artefacts and 



(a) Correct sampling. (b) Inappropriate sampling. 


Figure 3: The effects of inappropriately sampling the restor¬ 
ing beam on image quality. These images are of the same 
field, centred on 3C295, and they share the same physical 
and colour scales. At left, the image is correctly sampled, 
with 3 pixels across the minor axis of the restoring beam. 
To the right, only 1.3 pixels have been used: this image is 
unsuitable for pipeline processing. 


spurious source detections, as illustrated in Fig. [3] The 
TraP aims to be robust against this, regardless of the origin 
of the images. 

The shape of the beam is also considered. Beam shape 
is influenced by a variety of factors, including the array 
configuration, pointing direction, observing frequency and 
integration time. A measured shape which is significantly 
at variance with expectation is indicative of a failure in 
imaging or calibration. 

TraP users can predetermine image rejection thresholds 
for over sampled and highly elliptical restoring beams as 
these may be observation or telescope dependent. All 
images with restoring beams which are under sampled (< 2 
pixels across the FWHM) are automatically rejected by 
the TraP. 


4-2.3. Check for proximity to bright radio sources 

Poor subtraction of bright sources close to the target, 


during either ‘demixing’ (van der Tol et al. 2007) or subse¬ 


quent calibration can lead to residual structures or elevated 
noise levels in the resultant images. Problems are typically 
observed close to the Sun, Jupiter and the so-called “A- 


Team” of bright extragalactic radio sources (de Bruyn et al. 


2009): their extremely high fluxes may cause issues within 


target fields up to several tens of degrees away, depending 
on the observing configuration. To mitigate these effects, 
the TraP rejects images where the angular separation be¬ 
tween the target field and a bright source is less than a 
user-specified threshold. 

4-3. Source detection and measurement 

The TraP uses a custom-developed source detection 
and measurement system (“sourcefinder”). The algorithms 
implemented are partially based on those used in SExtrac- 


tor (Bertin and Arnouts 1996), but have been extensively 
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re-worked and extended, most notably to provide least- 
squares fitting of detected sources with elliptical Gaussians 
and a rigorous handling of the correlated noise properties of 
radio images. In addition, it provides an implementation of 


a false detection rate algorithm (Benjamini and Hochberg 


1995 Hopkins et al. 2002). 


In brief, given an input image, the system performs the 
following procedure: 


1. Model and subtract a variable background across the 
image; 

2. Model the RMS noise across the image; 

3. Identify islands of contiguous pixels which appear 
above some multiple of the RMS; 

4. Decompose multi-component islands into their con¬ 
stituent parts; 

5. Perform initial estimation of source parameters; 

6 . Fit each component with an elliptical Gaussian and 
return the result. 


This process results in a list of measurements describing 
the parameters of each source, including its position, flux 
density and shape as well as the significance of the detection: 
see Table [I] for details. 

More detail on each of the sourcefmder stages is given 
below. For a thorough treatment, the reader is referred to 


Spreeuw (2010), while Carbone et al. (in prep.) presents 


the results of extensive testing of the sourcefmder system. 


where F(na ) is the cumulative distribution function for 
the assumed normal distribution over the range [— oo,n<r]: 


F(na) = 0.5 + 0.5 x erf (2) 

= 0.5+ 0.5 x erf (3) 

2 f n 2 

= 0.5 + 0.5 x -— / e-* 2 d t. (4) 

V 77 Jo 

Inverting this, the threshold for clipping is 

na = aV 2 x erfc -1 ( ——- ] (5) 

\ 2iVindep / 


where erfc 1 is the complementary inverse error function 


the sample variance based on this indepen¬ 
dent pixel counlQ That is, 


(Gautschi 1972 


We estimate 


a 


2 

me as 


N inde P 
-^Mndep 1 ' 


( 6 ) 


where cr r 2 neas is the measured sample variance and x repre¬ 
sents individual pixel values. However, note that measuring 
the variance of a sample which has been clipped at some 
threshold T causes us to underestimate the variance as 
follows: 


4-3.1. Background and RMS estimation 

The background and noise characteristics of radio as¬ 
tronomy images are complex. In particular, noise across 
the image is correlated due to the nature of the imaging 
process. Since our source detection procedure relies on 
identifying pixels above the RMS noise, careful modelling 
of the background and noise maps is essential. 

We start by dividing the image into a rectangular grid. 
The dimensions of the grid cells are user-specified; they 
should be chosen such that they are significantly larger 
than extended sources visible in the image, but fine enough 
to trace background variations. 

We estimate the background and RMS in each grid cell 
through a process of iterative clipping around the median 
to reject source pixels. While doing this, it is important 
not to bias the result by rejecting source-free pixels, and 
to take account of the noise correlation scale. 

On each iteration, we label the total number of pixels in 
a given cell N. We then define the number of independent 
pixels as 7Vinde P = N/Nd ep , where N ^ ep is the number of 
pixels per synthesized beam. We assume that source-free 
independent pixels follow a normal distribution, while pix¬ 
els contaminated by sources do not. We therefore reject 
all pixels that fall more than some threshold na from the 
median, where a is the standard deviation of the distribu¬ 
tion. The value of n is chosen such that we will reject, on 
average, one half of one normally distributed pixel. That 
is, 

-/Vmdep x 2 x (1 - F(na)) = 0.5 (1) 


s 


^/ T T £ 2 exp(^)dx 


775^ f-T ex P ( 2 ^) dx 
We invert this to estimate a corrected variance 


2 2 
a — a.~. 


-\/27rerf (Tja\J 2) 


’ \f2neri (T / a^/2 ) — 2 T exp (T 2 /2a 2 ) /er 


( 7 ) 


( 8 ) 


Note that the ratio, T/a , of the clipping limit to the under¬ 
lying standard deviation is simply the value of n which was 
derived following Eq. [l]in the previous iteration. Finally, 
following Bolch (1968), we apply a further correction to 


estimate the standard deviation, a, as 


" = = <9) 

At this point, if any pixel values are more than the calcu¬ 
lated na from the mean, they are removed from further 
consideration and a new iteration is started. Otherwise, 
the clipping process has completed. 

After the pixel values have been clipped, if the field is 
crowded (|mean — median|/(std. dev.) > 0.3) we take the 


8 An alternative approach would be to calculate the pixel autocor¬ 
relation function and the corresponding estimation bias correction 
factor (W olter| , 1 984| ); the practical difference is minimal for plausible 
values of -/V;ndep- 
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background as equal to the median; otherwise, we estimate 
it as 


background = 2.5 x median — 1.5 x mean (10) 

following |Bertin and Arnouts 

Background and RMS maps are then derived by bilinear 
interpolation of the background and standard deviation 
calculated in each grid cell. 


4-3.2. Identification of significant pixel “islands” 

The background map is subtracted from the input image. 
The noise map is multiplied by a user-specified constant 
(“n”) to provide a threshold for finding sources. Sources 
are identified where the value of the background-subtracted 
input data is larger than the detection threshold map. 

Imposing an na threshold in this way implies that the 
fraction of pixels falsely labelled as sources in the image 
will be smaller than 



It may be more convenient, however, to quantify the 
number of false detections independently of the number of 
pixels processed. It is possible to control the false discov¬ 
ery rate using the algorithm described by |Benjamini and| 


Hochberg (1995). This permits the user to specify a maxi¬ 


mum number of “false positives” (noise peaks erroneously 
identified as sources), and the algorithm automatically 
chooses an appropriate threshold for identifying sources. 


4-3.3. Deblending 

After the peaks of pixel islands have been identified, 
the islands are expanded to consist of all pixels which are 
contiguous with the peak pixels and which are above a 
further user-specified constant (“a”, with a < n) times the 
threshold map. 

These expanded pixel islands may consist of a sin¬ 
gle source, but may also include two or more sources in 
close proximity. We separate out the components of these 
composite sources in a process referred to as “deblending” 


Bertin and Arnouts 

1996 

Spreeuw 

2010 

Hancock et al. 


2012). In the deblending process, we define a number of 


sub-thresholds exponentially spaced between the lowest and 
highest pixel values in the islands. We iterate through the 
thresholds from lowest to highest, at each threshold check¬ 
ing if the island has split into two (or more) non-contiguous 
components, each containing a significant fraction of the 
integrated flux density, and each with a peak above the 
detection threshold na. If these criteria are met, we split 
the island and regard it as two separate sources. Both the 
number of sub-thresholds and the fraction of the integrated 
flux density required for a sub-island to be regarded as 
significant are user-set parameters. 


4.3.4. Estimating source parameters 

The peak flux density, P max , can be approximated by 
the value of the maximum pixel in the island. However, 
the true source peak will not coincide with the centre of 
the pixel. Therefore, we extend the method described by 


Rengelink et al. (1997), based on the assumption that the 


true peak lies at a random position within the pixel. This 
results in a correction factor of 


M2) 






dxdy (12) 


where M, m and 6 are respectively the major and minor 
axes and the position angle of the synthesized beam and 
the integral runs over the pixel. This correction factor is 
multiplied by P max to produce the output peak flux density. 

The total flux density, P, is simply the sum of the pixel 
values 

f = ( 13 ) 

ies 

where is the value of the pixel at position Xi, yi and i £ S 
indicates all the pixels in a particular island. 

The position of the centre of the island is given in pixel 
coordinates as: 


x,y = 


SzgS Ii x i 


SigS liVi 
SigS F. 


(14) 


The position angle of the semi-major axis, measured counter¬ 
clockwise from the y-axis, is given by 


tan(20) = — Xy 
x 2 - y 2 


(15) 


The semi-major (M) and semi-minor (to) axis lengths are 
initially estimated as 



These axis lengths are underestimated due to the aa cut at 
the edge of the island. They are corrected by multiplying 
by a factor (1 + ln(T/P max )/(P max /T - l)) -0 ' 5 , where T 
is the value of the RMS map times the analysis threshold 
a at the pixel position of P max . 


4-3.5. Gaussian fitting 

An elliptical Gaussian is fitted to each island by min¬ 
imizing the error function using a modified Levenberg- 
Marquardt algorithm (More, 1977). By default, the es¬ 
timated source parameters calculated above are used as 
initial values for fitting, and all parameters are allowed to 
vary. However, the user may optionally choose to hold one 
or more parameters to a fixed, user-specified value. Typi¬ 
cally, this is used to constrain the fitted shape of point-like 
sources to that of the restoring beam. 

Uncertainties on the parameters xo,yo (the fitted pixel 
position), a, S (the position in celestial coordinates), 9m , 9 m 
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(the lengths of the major and minor fitted axes), C (the 
fitted peak flux density), I (the integrated flux density) 
and 0 (the position angle) are calculated following Condon 
(1997), Condon et al. ( 1998| ) and Hopkins et al. (2003). We 
start by defining a generalized “signal-to-noise ratio” as 


P 2 = 



\ 1+ °b] 

otM 

[i+ 0fe2 l 

A9b 9 b 

[ 9 2 m \ 


Id 


r 2 

CT (17) 


where a is the RMS noise at the location of the source 
and aM,Oi m = (1-5,1.5) for amplitude errors, (2.5, 0.5) for 
errors on x and 9 m and (0.5, 2.5) for errors on y, 9 m and 


(j) (Condon 1997). 


Given the above definitions, we apply the relationships 


described by Condon (1997) and Hopkins et al. (2003) to 
obtain 


C 2 

= 81 n 2 pF 


M 


CTa 

VM 

4 (O' 


= 81n2-^ 

u m 

°~L 

° 2 m 

- 0 2 m? 


2 

2 




= °l 0 sin2 < 

= < COS 2 ( 


°i 0 cos2 - 

- a 2 yo sin 2 « 


(18) 

(19) 

( 20 ) 

( 21 ) 

( 22 ) 

(23) 


Note that, for simplicity, the above assumes that the fitted 
major axis of the source aligns with the y axis of the pixel 
grid. If this is not the case, a further rotation is required. 
Finally, the integrated flux density and its variance are 
given by 


I = C 


9 M 9 m 

SB^b 


a_i _°c_ , 9 b 9b 
P~C 2 9 M e m 


0 2 m 


Z_Brn 

9 2 


(24) 

(25) 


Name 

Units 

Notes 

Right ascension (a) 

J2000 deg. 

Declination (6) 

J2000 deg. 

Peak flux density 

Jy/beam 


Integrated flux density 

Jy 


Significance 

- 

Peak/RMS 

Lengths of semi-axes 

arcsec 


Position angle 

deg. 

North through east 

Error radius 

arcsec 

Absolute error on 



source centroid 


Table 1: Parameters returned by the sourcefinder routines. 


4-4- Source association 

Each individual astronomical source detected in a given 
set of images is assigned a unique identification in the 
form of an entry in the “running catalogue”. The running 
catalogue ties together a series of measurements made in 
individual images with aggregated information about the 
source derived from those measurements (its position, vari¬ 
ability information, etc; see 14.7). The complete set of 
measurements associated with a particular running cata¬ 
logue entry comprise its lightcurve. 


4-4-1- Association procedure 


The association procedure adopted is based on|de Ruiter 


et al. 

1977 

, Sutherland and Saunders 

(1992 

and Rutledge 

et al. 

(2000 

), and is described in detail in Scheers 

(2011). 


For each measurement, the source association procedure 
searches for counterparts in the running catalogue. The al¬ 
gorithm relies on the de Ruiter radius, the angular distance 
on the sky between source i and its potential association 
counterpart j normalized by the positional error of both 
sources. The de Ruiter radius is defined as 


' »>3 



ay ) 2 cos 2 ((Sj + Sj)/2) {6j - 6j) 2 

a li + a l + 


(31) 


After fitting, the Gaussian restoring beam is decon¬ 
volved from the resultant source parameters using an algo¬ 


rithm derived from that provided by AIPS (Greisen 2003) 


so that the deconvolved shape parameters 9m, 9 n 
are given by 

P 2 = (&M ~ Qm) + (&B - 9b) 

- 2(6 2 m = 9 2 J(9 2 b -9l) cos 2(0-$) 

9m) ~ (9% + 9b) + P 

= ( 9m + 9m) — (9% + 9 2 ) — P 


2 9 2 m — [9\ 


1 M 


1 


p = -atan 




(^-e)cos2(<),-$)-(0|-0 2 ) 


and ip 

(26) 

(27) 

(28) 

(29) 

(30) 


Following this procedure, the parameters described in 
Table [T] are returned by the sourcefinder and are ready for 
insertion into the pipeline database. 


where a n is the right ascension of source n, S n is its decli¬ 
nation, and o q represents the error on the quantity q. 

If sources i and j are genuinely associated, their posi¬ 
tional differences will be due to measurement errors, and 


hence follow a Rayleigh distribution (e.g. de Ruiter et al. 


1977). The probability of source association at r > p is 


then 


pOO 

Pr{r >p)= r exp(—r 2 /2)dr = exp(— p 2 /2). 
J r—p 


(32) 


This may be used for determining the search radius, r s , 
of the area that will be scanned for possible counterparts: 
a search radius of r s < 3.71, will miss a factor of ICC 3 of 
the possible association counterparts, while r s < 5.68 will 
miss a factor of ICC 7 . 

Given the above definition, the source association proce¬ 
dure regards a particular measurement as being associated 
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with a given running catalogue source if their positions are 
no further apart than the semi-major axis of the restoring 
beam and the de Ruiter radius is less than a user-specified 
threshold. Note that the calculations above only consider 
repeat measurements of a single, isolated source. If the 
TraP is to be used in processing observations of crowded or 
extremely transient-rich fields of view, this will require fur¬ 
ther consideration of the trade-off in search radius between 
missed self-associations, and spurious associations between 
distinct sources. Making an optimal choice of search radius 
for a crowded field will depend on the precise spatial clus¬ 
tering of sources, an issue which is not investigated further 
here. 

4.4.2. Association types 

The procedure described above does not guarantee a 
unique association between each source measurement and a 
single running catalogue source. Instead, there are various 
possible association topologies, as illustrated in Fig. [4] and 
discussed below. 

Note that the association is based upon the positions 
and associated uncertainties of a particular measurement 
and the equivalent aggregated quantities for running cat¬ 
alogue sources; no reference is made to time or frequency 
in assessing whether a particular pair is associated. The 
discussion below refers to “time”, but the same consider¬ 
ations apply to association along a frequency axis. This 
has the consequence that the order in which data is asso¬ 
ciated affects the result, and hence the reproducibility of 
a particular analysis. This is discussed in more detail in 

fro 

No association. The source measurement cannot be asso¬ 
ciated with any previously catalogued sources. We regard 
this as a newly detected source, and create a catalogue 
entry for it. 

One-to-one. The flux density measurement is unambigu¬ 
ously associated with a single running catalogue entry, as 
shown in Fig. [4a[ The average flux density of the source 
Li is / 1 ... 4 . 

Many-to-one. Many-to-one associations exist when two or 
more running catalogue sources satisfy the association crite¬ 
ria for a given measurement. We record this in the database 
as multiple separate one-to-one associations. This is illus¬ 
trated in Fig. [4b} at times t-\ and t 2 distinct lightcurves, L\ 
and L 2 are being tracked. However, at 1 3 a single source 
is detected which can be associated with both of these 
lightcurves. This could happen, for example, if the third 
observation was made at a lower angular resolution. 

Note that the single flux density measurements and 
fe are now independently included in two separate sources: 
L 1 has average flux density / 1 , 3 , 5,6 and L 2 has average flux 
density f 2 ,A, 5,6 ■ The total brightness reached by summing 
all catalogued source flux densities has been artificially 
increased. Since all individual measurements are stored, 



(c) One-to-many. 



(d) Many-to-many. 



(e) Reduced many-to-many. 

Figure 4: Types of source association, corresponding to 
those described in 14.4.2 Flux density measurements taken 
at a particular position at time t a are labelled fb- The asso¬ 
ciation procedure knits flux density measurements together 
between timesteps to form lightcurves which are identified 
with particular running catalogue entries identified as L c . 
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it is possible for users to correct for this in their analysis 
according to their particular requirements. Future versions 
of the TraP may support the division of flux density from 
a single measurement among multiple sources (j|10|). 


One-to-many. In the reverse of the previous case, a sin¬ 
gle running catalogue source can be associated with more 
than one source measurement from a given image. This is 
handled by splitting the catalogue source into two indepen¬ 
dent sources, L\ with average flux density /i, 2 , 3,5 and L 2 


with average flux density /i, 2 , 4,6 as shown in Fig. 


4c 


As 


in the many-to-one case, some source measurements are 
included in multiple lightcurves, artificially increasing the 
total brightness of the catalogue. 


Many-to-many. A many-to-many association occurs when 
more than one running catalogue source can be associated 
with more than one extracted source. If the procedures 
described above were applied, every possible combination 
of catalogue sources and new measurements would result in 
a new lightcurve: the database complexity would increase 
quadratically, and the situation rapidly becomes untenable. 
This is illustrated in Fig. |4d[ To prevent this, many-to- 
many associations are reduced to one-to-one or one-to- 
many associations by choosing only the source pairs with 
the smallest de R.uiter radii. By applying this procedure, 
the association topology is reduced to a simpler case such 
as that shown in Fig. [4e] rather than the eight separate 
lightcurves produced by “pure” many-to-many association, 
we are left with just three. 


4-5. Null-detection handling 

We use the term “null-detection” to describe a source 
which was expected to be observed in some particular 
image—it had been observed in previous images of the 
same field with a flux density above our detection threshold 
for this image—but which was not detected by the source 
finding procedure (f 4.3). 


After the association step (14.4), we build a list of all 
sources which: 


• Exist in the list running catalogue of known sources, 
having been observed (in any frequency band) at an 
earlier timestep; 


catalogue entry; the measurement is marked as being due 
to the null detection procedure. 

In the current release of the TraP, once a running cata¬ 
logue entry has been created, forced fitting at that location 
will always be attempted in the absence of a matched 
blind-detection. Ultimately, an accumulation of marginal 
single-epoch false detections due to noise will cause a very 
large number of irrelevant forced fits to be performed and 
stored. This may be mitigated by expiring these marginal 
detections from the list of positions to be measured if they 
are not re-observed after a certain period of time or num¬ 
ber of revisits. Automatic expiry according to user-defined 
criteria will be included in a later release of TraP. 


4-6. Monitoring list sources 

We anticipate that there may be some sources which 
are important to particular science cases which may not 
always be detected by the default sourcefinding procedures 
employed by the TraP. It is therefore possible for the end- 
user to supply a “monitoring list” of positions at which a 
measurement will always be made, regardless of the results 
of the sourcefinding step. The current version of the TraP 
assumes that input images have been correctly registered 
when making these measurements: it makes no attempt to 
correct for astrometric shift. 

For each image which covers the location of a position on 
the monitoring list, a measurement is taken at its location. 
The same procedure is used as for null detections ((4.5): a 
modified version of the algorithms described in §4.3| which 
hold the position of the measurement constant. 

For each monitored position, a running catalogue source 
is generated which contains only a chain of one-to-one 
associations of measurements at the user-specified position. 
Sources monitored in this way are not included in the 


general association procedure described in (4.4.1 


4-7. Aggregate source properties 

For each entry in the running catalogue, we now have at 
least one individual measurement (corresponding to the cur¬ 
rent image) and potentially others corresponding to other 
images (presumably representing observations at different 
times or frequencies). We combine these measurements to 
estimate the true source properties as follows. 


• Were not associated with any source detected in the 
current image. 

For each of these null detections, we use a modification 
of the sourcefinding procedure: the same techniques are 
employed for background and noise estimation and source 
measurement as previously described, but, rather than 
being based on islands of pixels above some threshold, 
the peak of measurement is constrained to be fixed to 
the catalogued position of the null detection. No other 
additional constraints are placed on the measurement. 

After the measurement has been made, it is stored 
as a one-to-one association with the appropriate running 


4-7.1. Mean values 

For each property of a given source we store both the 
arithmetic mean and a weighted mean. For a series of 
measurements of property x , we denote the arithmetic 
mean of x as x. We define the weight of a measurement of 
property x as 

w x = 1/al (33) 

where a x is the uncertainty on that measurement. The 
weighted mean of N such measurements is then 


4 * 


lsi=l w XjXj 

\~~\N 

Ei=1 w *i 


(34) 


12 










Using these definitions, for each source we calculate 
both the arithmetic and weighted means of: 

• The source right ascension (J2000 degrees); 

• The source declination (J2000 degrees). 

For each frequency band for each source, we calculate 
both the arithmetic and weighted means of: 

• The peak flux density of the source in this band 
(Jy/beam); 

• The integrated flux density of the source in this band 

(Jy)- 


4-7.2. Variability metrics 

After source association, each running catalogue entry 
corresponds to a multi-frequency lightcurve of a particular 
source. We search for transient and variable sources by 
examining the properties of these lightcurves. 

For each source, we label a measurement of its flux 
density at a particular frequency, as I v with uncertainty 
ay. Based on these quantities, and using the same notation 
as above, we define the flux density coefficient of variatiorj^] 
over N measurements as the ratio of the mean flux density 
to the sample standard deviation s, thus: 



Using the same definition of weight w v = 1 /c 2 as above, 
we can also express the significance of the flux density vari¬ 
ability. We make use of the familiar reduced-y 2 expression 
in conjunction with the weighted mean flux density, 


Vv = 


1 

N ~ l k <i 


(36) 


For a given ? 7 „, the probability that we can reject the 
null hypothesis—that the source under consideration is not 
variable—is given by 


variable — 1 


p{jiv',N - 1) dr]/. (37) 


where p(x, n) is the y 2 probability density function for x 


over n degrees of freedom (see, for example, Kesteven et al. 


1977). 


V v and ip, are calculated and stored for every lightcurve 
whenever a new source measurement is added. Since vari¬ 
ability metrics are stored per association, we can track how 


“This quantity is occasionally referred to as the ‘modulation index’ 


in astronomical literature (e.g.|Narayan 

1992||Gaensler and Hunstead 

2000[ | Jenet and Gil 2003 Bell et al. 

2014|. The present authors 

prefer ‘coefficient of variation’ due to its more widespread use ( McKay 

1932; Hendricks and Robey, 1936| Lande||l977| fFreund and Williams 


2010 ) and because it avoids any possible confusion with other fields 
(e.g. IWhitaker 19961. 


the variability parameters of a source have changed with 
time. This is particularly useful in the case of those sources 
which have shown evolution in their behaviour over time. 

It is worth noting that it would be relatively straight¬ 
forward to extend the TraP to support the calculation 
and storage of other variability metrics beyond the two 
described above. It is expected that extended testing and 
experience in processing data from various sources will 
guide future development in this area. 

4-7.3. Reproducibility of results 

Reproducibility of pipeline results is of paramount im¬ 
portance: the end user should be confident that repeatedly 
running the TraP on a given dataset with the same config¬ 
uration should always produce the same results. This has 
important consequences for the association and aggregation 
procedures. For example, consider a particular running 
catalogue source R and two source measurements, M\ and 
M 2 , taken from different images. If Mj is inserted first, 
it is associated with R. On association, a new aggregate 
position for R is calculated (which may or may not be 
consistent with association with M 2 ). On the other hand, 
if M 2 is inserted first, the resulting aggregate position for 
R is not consistent with association with M\. In short, 
the order in which the images are processed influences the 
association procedure, and hence changes the outputs. 

In order to mitigate this effect, the TraP only guar¬ 
antees reproducibility of output if the input is in mono- 
tonically increasing order of time. If two or more images 
with different frequency or Stokes parameters but the same 
timestamp are processed, the TraP will automatically sort 
them along these axes before processing. This is not, in 
general, possible along the time axis, which is potentially 
unbounded. 


5. Data products 

After all the stages described in (JdJhave been completed 
for a given image cube, the core pipeline loop is complete. 
The complete running catalogue, together with ancillary 
data describing the images which have been processed, 
pipeline configuration and other metadata, is now stored 
in a database, the structure of which is described in detail 
in §6.4[ At this point, the pipeline may start processing 
the next image cube. Simultaneously, the contents of that 
database may be used to support scientific analysis and 
alert generation. 


5.1. Identifying transient events 

As described in (|3j we distinguish between newly de¬ 
tected and variable sources. Both are scientifically signifi¬ 
cant and may form the basis for a “new transient” alert. 
Depending on context and science goal, these alerts may 
simply result in a particular source being brought to the at¬ 
tention of the pipeline end user, or they may be distributed 
more widely to the community. The technical mechanism 


used for alert dissemination is described in 16.6 
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5.1.1. New detections 

As described in §4.4.2| a measurement which cannot be 
associated with a previously catalogued source is regarded 
as a new detection. Such a new detection might correspond 
to a new transient source appearing, but it could also 
simply mean that this area of sky had not previously been 
surveyed to a depth adequate to detect this source. 

In order to distinguish between these possibilities, the 
TraP keeps track of the fields-of-view and sensitivities of 
all images it processes. When a new source is detected, its 
flux density is compared against the recorded sensitivities 
of images covering the same area to see if it could have 
been detected previously. If so, it may be regarded as a 
new transient. 

In practice, and as described in 14.3.1 noise (and hence 
sensitivity) is not constant across any given image. It is 
possible that a particular source could have been detected 
if it fell in a low-noise area of a previous image, but not if 
it fell in a high-noise area. We therefore record both the 
highest and lowest RMS values recorded in each previous 
image (<7 max ,i ,v and cr m i n> i tI/ in image i at frequency v) as 
well as the detection threshold (rq j!/ ) used when processing 
that image. For a flux density measurement we regard 
it as a marginal transient candidate if 


+ M) > I iiV > ^min ,j,i/ X (n jt „ + M). (38) 

where M is some user specified margin, applied to prevent 
stable sources with flux densities approximately equal to the 
detection threshold from being misidentified as transients. 
This marginal category would include both genuine but 
faint transients, and steady state sources which change in 
significance as noise varies between images. In the case 
that 

li.v ^ ^maxj> X ( TLj )Iy T A/) Vj <t i (39) 
we regard the source as likely to be a new transient. 


5.1.2. Variable sources 

As per §4.7.2 three separate variability metrics—the 
flux density coefficient of variation V v , the significance 
of the variability rj and the probability of variability 
-^variable—are stored whenever a new source association is 
made. We can therefore search the catalogue for sources 
which meet the following criteria: 


• V v is above a user-specified threshold; 

• rjjy is above a user-specified threshold; 

• ^variable is above a user-specihed threshold; 

• The number of source measurements in the lightcurve 
is above a user-specified threshold. 


Choosing appropriate values for these thresholds is 
a matter of user configuration and will depend on the 
details of the science case as well as the particular dataset 
under investigation. Section [8] gives an overview of possible 


considerations, while Rowlinson et al. (in prep.) presents 
a detailed consideration of how optimal parameters might 
be chosen. 


5.2. Database interface 

While automatic routines may be used to scan the 
database for transients as it is constructed using the meth¬ 
ods described in §5.1[ it is likely that many scientific dis¬ 
coveries will be based on expert analysis of the database 
contents. Here, we describe the systems by which this is 
made available to end users. 


5.2.1. Direct database queries 

The core database engine is one of various off-the-shelf 
data base management systems, as discussed in detail in 
{6A Given appropriate permissions^] the database may be 
queried directly using SQIp~| Expert users may write SQL 
scripts or use the command line interface to the database 
system to manually construct complex queries selecting 
exactly the information they require. 


5.2.2. Web interface 

While the ability to query the database for arbitrary 
information which answers a specific science question is 
undeniably powerful, it requires a significant investment 
of time on the part of the end user to become fluent in 
the relevant technologies. Further, it is often convenient 
to have an at-a-glance summary of the contents of the 
database and the important results. For this reason, the 
TKSP has developed Banana, a web-based interface to the 
database. Banana enables the user to conveniently examine 
the contents of the database, viewing details (including cut¬ 
out images) of all source measurements, plotting lightcurves 
of all sources in the running catalogue, selecting potential 
transients based on their variability metrics, and so on. 

Banana is open-source (it is released under a Modi¬ 
fied BSD license ^j) and is developed and released by the 
TKSP independently of the TraP. It is freely available for 
download*^! 


5.2.3. High volume archival data-mining 

For modest input dataset sizes (thousands of images, 
tens or hundreds of sources per image), the total volume of 
data stored is modest: on the order of, perhaps, gigabytes. 
However, as per §2.2| the TraP ultimately aims to support 
long term operation of the LOFAR RSM, which will be 
capable of producing thousands of images per hour. It is 
also a stated aim of the project to make the lightcurve 


10 Database permissions are controlled by the administrators of 
a particular TraP installation; it is possible for them to grant per¬ 
missions to both query and modify the data to arbitrary users if 
required. 

11 Structured Query Language, the de-facto standard for interacting 
with relational database systems. 

12 http://opensource.org/licenses/BSD-3-Clause 
ld https://github.com/transientskp/banana 
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archive available to the community as a legacy resource. 
Efforts are currently underway to both develop database 


systems capable of handling this volume of data (Scheers 


2011 Scheers et al. in prep.), and a batch query system 


akin, for example, to CasJobs (O’Mullane et al. 2005) is 
under consideration. Ultimately, we also hope to make 
data available using a Virtual Observatory-compliant in¬ 
terface. However, the currently-available IVOA model for 
time series data (Graham et al. 2014) is an intentionally 
minimal interim solution; we prefer to wait for more mature 
standards to become available (e.g. McDowell et al. 2015) 
before proceeding. 


6. Implementation 

In this section, we turn our attention to the underlying 
technology which implements the workflow described in 
sections |3] and 0] 


6 .1. Architecture 


The TraP is structured as a series of pipeline “steps”, 
each of which performs a logically discrete operation in 
turn. These steps, and the relationship between them, are 
shown graphically in Fig. [2] 

The operations carried out by the TraP naturally split 
into those which involve directly interacting with image 
data (such as loading data from disk or finding and mea¬ 
suring sources within it) and those which operate on mea¬ 
surements derived from the data (source association or 
aggregate calculation). 

Interacting with images and performing measurements 
upon them is most effectively accomplished by bespoke soft¬ 
ware packages developed by members of the TKSP which 
directly encode the required algorithms. We have developed 
a series of such packages using the Python programming 
language^ The choice of Python for this task, together 
with a discussion of the approach taken, is motivated in 

do 


The derived data which results from source measure¬ 
ments is highly structured. It can be efficiently stored using 
a relational database system (Codcl 1970). The contribu¬ 
tion of the database goes beyond mere storage, however: 
by performing calculations within the database itself it 
possible to operate on the entire corpus of data efficiently 
and with minimal overhead due to data transport. In this 
way, the database becomes the core computational engine 
of the TraP. The design and structure of the database is 
described in §6.4| 

A control layer sits above the core scientific logic as 
defined in Python code and the database. This control layer 
defines the structure of the pipeline—effectively connecting 
the components together in the correct order—as well as 
providing utility services such as parallelization and task 
distribution, which we describe in 96.31 


14 http://www.python.org/ 


Finally, and in addition to the pipeline routines de¬ 
scribed, the TraP offers the option to save a copy of all 
the pixel data processed to a separate document-oriented 
database for later use by the Banana web interface ( Q5.2.2 1. 
This is described in §6.5| 


6.2. Python 

Python is the primary programming language used in 
the TraP. We consider Python to be the default choice 
for astronomical software development where performance 
is not the critical consideration in the near to intermedi¬ 
ate future. It provides a flexible and expressive language 
together with a wide ecosystem of scientific and other li¬ 
braries, and it is easily extensible using code written in 
lower-level languages where maximum performance is re- 
quirecf^l Furthermore, thanks to projects like IPython 


(Perez and Granger 


2007) and AstroPy (Astropy Collab¬ 


oration et al. 2013), Python is also increasingly finding 


a role in the daily workflow of many astronomers as an 
interactive data analysis toolbox. Although we do not 
directly use these tools in the TraP, this familiarity then 
lowers the barrier to entry on larger projects as the novice 
coder becomes more proficient, potentially widening the 
pool of future maintainers and contributors to open-sourced 
scientific codes. 

Although we have had great success using Python, a sig¬ 
nificant downside is that, as a dynamically typed language, 
there is a risk of run-time type errors. We have countered 
this by adopting a strongly test-focused development style 
and building an extensive suite of TraP unit tests. This is 
a topic we return to in §7.4| 


6.3. Parallelization and distribution 

Some operations which are carried out by the TraP 
can be performed concurrently on multiple datasets. For 
example, the initial source finding and measurement step 
( |4.3[ ) can be performed on many independent images si¬ 
multaneously without changing its results. We can exploit 
this intrinsic parallelism to obtain the best possible run¬ 
time performance by distributing processing across multiple 
CPU cores (or even distinct machines) and scheduling as 
many operations to run concurrently as is possible. 

Of course, the performance improvement which may be 
achieved in this way is limited: some operations (such as 


source association, 14.41 cannot run concurrently (we rely 


on data being processed in a particular order to ensure 
reproducibility; j |4.7.3 ); there is an intrinsic ordering of 
pipeline steps (it is impossible to perform source association 
before source measurement is complete); and, for reasons 
of reproducibility, we mandate that all steps relating to 
images corresponding to a particular observation time are 
complete before a subsequent observation time can start 


15 We have experimented with writing portions of the TraP in 
Cython (http://www.cython.org/) for performance reasons with 
some success, but no Cython code is shipped with the current release. 
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processing (f 4.7.3). These intrinsically sequential parts 


of the processing limit its overall performance (Amdahl 


1967). 


We have implemented the TraP in such a way that the 
definition of the underlying algorithms is independent of the 
method used to start tasks and collect results. In this way, 
it is possible to insert different task scheduling back-ends 
which support different parallelization and distribution 
techniques. Three are currently supported by the TraP: 


• The serial back-end runs tasks sequentially using a 
single Python interpreter. Using the standard Python 
interpretei]^] means that all (non-database) process¬ 
ing takes place in series on a single CPU core. 

• The multiproc back-end uses the multiprocessing 
packag(pl to schedule jobs on multiple CPUs within 
a single machine concurrently. 

• The celery back end uses CelerJ^] an asynchronous 
task distribution system, to marshal the distribution 
of concurrent TraP tasks across a cluster of multiple 
machines. 


The end user may select which back-end to use when 
invoking the TraP from the command line. 

Note that the celery system does not arrange for data 
to be transmitted across the cluster. If, for example, it 
is used to distribute a source finding step across multiple 
images, it is required that each machine have access to the 
particular images which it is to process (perhaps on its local 
disk or on shared storage). This is a convenient match to 
the imaging process on the LOFAR cluster, which deposits 
image data on the cluster node which was responsible for 
creating it. We have, nevertheless, prototyped an image 
transmission system which is better integrated with the 
TraP, but it is not included in the current release. 

Many-core exploitation works well for the ‘embarrass¬ 
ingly parallel’ problem of sourcefinding across many dif¬ 
ferent images (corresponding to different frequencies and 
pointings). However, the bottleneck for processing a timestep 
then becomes the database operations, placing stringent 
performance requirements on the combination of query 
complexity and database back-end used, as covered in Sec¬ 
tion El 


the TraP may write to a given database during processing 
(but the database may be accessed by arbitrary read-only 
queries). However, both of the supported database manage¬ 
ment systems, as described in §6.4.1 provide for multiple 
isolated databases being hosted within a single system, so 
supporting multiple pipeline instances is straightforward. 
This could be used, for example, to support multiple inde¬ 
pendent monitoring campaigns. 

6 .4-1- Database management systems 

The TraP is developed and tested using two relational 
database management systems (RDBMS): MonetDEp^l and 
PostgreSQIp 7 ] 

MonetDB is an ongoing project to build a database with 
exceptional performance and scalability based on research 


into data organization and query optimization (Idreos et al. 


2012). The potential performance benefits of MonetDB are 


impressive, particularly when considering the ultimate data 


volumes expected from the RSM ((2.2) and, later, from the 


SKA. However, its development is driven by fundamental 
database research and scientific user groups, and it may 
occasionally perform in unexpected or undesirable ways. 
We therefore also verify the correct operation of the TraP 
and, where necessary, provide “production grade” deploy¬ 
ments using PostgreSQL, which has a long pedigree as an 
industry-standard database management system. 

The TraP pipeline code and the Banana web interface 
send queries to the database using SQL. Although SQL 
is a standardized language^] there is significant variation 
in its implementation between different database vendors: 
code that is written and tested against one database may 
unexpectedly fail when run on another system. Therefore, 
while the code in the TraP is designed to be standards 
compliant and database vendor agnostic, it is occasionally 
necessary to add special cases to work around different 
SQL dialects. To accommodate this the TraP provides a 
simple templating system for SQL queries. For example, 
we can accommodate both the PostgreSQL and MonetDB 
syntaxes for defining a function within the database as 
follows: 

{°/o ifdb monetdb 7,} 

CREATE PROCEDURE BuildFrequencyBands() 

{7. endifdb 7} 


6-4- Database 

All source measurements, together with metadata de¬ 
scribing the images from which they were taken, are stored 
and processed in a relational database. The contents of the 
database is itself one of the core products of the TraP (f 3.21, 
while queries run over the database are an intrinsic part 
of regular pipeline processing. Only a single instance of 


{7. ifdb postgresql 7o} 

CREATE OR REPLACE FUNCTION BuildFrequencyBands() 
RETURNS void 
AS $$ 

{7 endifdb 7} 

BEGIN 

— Definition elided 
END; 


1(, The reference “CPython” implementation. 

11 multiprocessing is part of the Python standard library. 
lf Tttp: //mu. celeryproject. org/ 


19 http://www.monetdb.org/ 

21 http: //www.postgresql. org/ 
21 iSO/itXJ 9075-1:2011 
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In this way it is easy to extend the TraP’s database support 
to encompass other systems if required. 

6-4-2. Database structure 

The structure of the current version of the TraP database 
is elaborate, consisting of many separate tables, some with 
tens of individual columns, complex relationships between 
them, and a variety of stored procedures. It is described 
in detail in the TraP 1 himlboolj^j A simplified version is 
shown in Fig. [5] 

It is expected that, as TraP development progresses, the 
database will evolve to meet new requirements. For this 
reason, the database is versioned, and the TraP does not 
support mixing results across versions. The current code- 
base does not offer explicit support for “schema migrations” 
(converting a database from one version to another without 
losing data), although versions with this functionality have 
been prototyped, and off-the-shelf schema migration toolj^] 
may be applicable. 

It is worth noting, for performance reasons, the database 


is not fully normalized (Kent 1983). 


6-4-3. Positional database queries 

Perhaps the most common operation required of the 


database, both during pipeline processing (e.g. (4.4) and 
during later analysis is the “cone search”: finding all ob¬ 
jects of a particular type (such as source measurements, 
running catalogue entries) within a given search radius of 
a particular position. Since this operation is so common, 
we give it special consideration. 

Given a particular target (a T , S T ), search radius r and 
list of positions (op, 6i), («2, 62)1 Sn), the simplest 

approach to finding all the positions which fall within the 
radius of the target is to iterate down the list calculating 
for each one the great-circle distance 

r T n = arccos(sin S T sin 6 n + cos S T cos S n cos \a T — a n \) 

( 4 °) 

and selecting only those sources for which r T n < r. While 
conceptually simple, this involves multiple calculations for 
every source to be checked: it is prohibitively computation¬ 
ally expensive given a large source list. 

Since this simple approach is impractical, we adopt an 


approach based on that described by Gray et al. (2006 


24 


This involves a hierarchical approach, filtering the list of 
candidate sources first on declination then on right ascen¬ 
sion before selecting candidates based on a Cartesian dot 
product. 

First, when any new source measurement is inserted 
into the database, or when the weighted mean position 


“-http://docs.transientskp.org 

23 For example, Alembic http://alembic.readthedocs.org/ 
24 Note that a variety of alternative approaches were considered, 
such as HTM jSzala.y et al.|[2006^ and Q3C l |Koposov and Bartunovj 
2006); practical experience and compatibility with the database man¬ 
agement systems informed the approach taken. 


((4.7.1) of a catalogue source is updated, we calculate a 
corresponding position on the unit sphere: 


x = cos S cos a 
y = cos S sin a 
z = sin 6 


(41) 

(42) 

(43) 


At the same time, we define a function zone(<5) which 
maps a declination to a particular “zone”, corresponding 
to a strip of the sky. Zones must increase monotonically 
with declination: zone(S m ) includes all declinations falling 
between zone(<5 m _i) and zone(i5 m+1 ). In the current ver¬ 
sion of the TraP we use the largest integer smaller than 6 
as the zone, thus: 


*(«) = !*]• 


(44) 


However, this definition is flexible: given the constraint 
above, future versions could adopt a zone definition with a 
more fine-grained resolution or variable zone heights (e.g. 
chosen to provide zones of uniform area). 

The Cartesian coordinates and the zone are stored in 
the database and are henceforth available for each running 
catalogue source and source measurement with no further 
run-time calculation. 

It is next necessary to describe the ‘inflation’ of angular 
distances in right ascension with declination. For example, 
at a declination of 0°, a circular region with radius 6 centred 
on right ascension a includes RAs in the range [a — 6, a + 9\, 
whereas at a declination of 90° it covers the complete circle. 


Following Gray et al. we define the function alpha(0, S) as 


alpha(0, 6 ) = 


arctan 

180 


sin 6 


yj cos(<5— 0) cos(<5+0) 


if \6\ +6 < 90; 
otherwise. 


( 45 ) 

In general, such a circle at arbitrary a, S can be said to cover 
the range [a — alpha)!?, 6 ), a + alpha(0, <5)]. alpha(0,6 ) is 
implemented as a stored procedure directly in the database 
so that it can be calculated for arbitrary 6 and <5 with 
minimal overhead. 

These definitions made, we are now able to quickly 
filter the list positions to be searched. First, we cal¬ 
culate the maximum and minimum zones in which tar¬ 
gets may be found, rejecting all those targets for which 
zone((5 n ) lies outside the range [zone(5 r _ r ), zone(5 T+r )]. 
Then we reject all targets for which a n lies outside the 
range [a T — alpha(r, S), a T + alpha(r, <5)]. By ensuring that 
the database is appropriately indexed on zone(<5„) and a 
this filtering can be done extremely fast. 

The above filtering reduces the potentially large list 
of positions to match to a much more manageable size. 
For each position, we now check whether it lies within 
the required angular distance of the target. Rather than 


calculating the great circle distance (Eq. 40) it is more 


efficient to use a scalar product based on the Cartesian 
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Figure 5: Simplified version of the TraP database structure. Each box represents a database table, while lines represent 
key-based cross-references. A single line refers to a unique key, while a split line indicates that a many-to-one relationship 
is possible. For example, an Image is a member of a single Dataset, but a Dataset may contain many Images. 


positions calculated according to Eqs. piT| pTTT| Thus we 
check for: 


x T ■ x n + y T ■ y n + z T ■ z n < cos r (46) 

Since (x n , y n , z n ) for each candidate is already stored in 
the database, the total amount of computation (and hence 
running time) is kept to a minimum. 

The above procedure fails in the case of a discontinuity 
in RA: at the meridian, we jump from 359° to 0°, breaking 
the check for a n lying within a given range. To work around 
this, if a source association query crosses the meridian, we 
rotate the RAs of the relevant sources by 180° to avoid the 
discontinuity, perform the association as normal, and then 
rotate the results back to the original orientation. 

6.4-4- Iteratively updating aggregate quantities 

As per §4.7[ we store weighted mean positions, flux den¬ 
sities and variability indices for all sources in the database. 
When a new measurement is appended to a particular run¬ 
ning catalogue entry, it would be possible to re-calculate 
these quantities from scratch by averaging over all the ex¬ 
isting information about that source in the database. This 
would clearly be inefficient, though. Instead, we update 
these quantities iteratively. 

The arithmetic mean of some property x over N mea¬ 
surements is 

1 N 

Xn = n Xi ( 47 ) 

V 2=1 

where Xi is the i th measurement of x. When a further 
measurement of x is taken, updating the mean iteratively 


is straightforward: 


Xn+1 = 


N x n + x N+1 
N + l 


(48) 


As per (4.7.1 we also calculate the weighted mean 
defined as in Eq. 34 The weight of measurement Xi is w Xi ; 


the sum of all weight over N measurements is 


N 


= E 


Wx 


(49) 


i=1 


Given £ XJV and W XN it is then possible to express the 
weighted mean after N+l measurements as 


£ 


xn+i 


W XN ^ XN T w Xn+ 1 x n+1 

W XN + w XN+1 

Nu>x n £, xn “ 1 “ w xn + 1- X N+1 

Nw Xn + w XN+1 


(50) 

(51) 


Therefore, by storing the number of measurements (N) 
and the average weight (w Xn ) in addition to the weighted 
mean, we can iteratively update the mean and weighted 
mean source properties as new measurements are added to 
the database without revisiting all previous measurements. 

As measurements are associated with catalogue sources, 
we keep track not only of the mean parameters but also of 
the evolving variability parameters, V v and 7y„, as described 
in ' jT7A| 

Updating V v is straightforward. We store both the 
mean and the mean square flux density per band, /„ and 
11 and update them using the procedure described above. 
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They can then be directly used to calculate V v based on 
Eq. [351 


To handle rj„, we substitute Eq. 33 into Eq. 36 to get 

N — 1 \ N 



N / 1 


N — 1 \ N 




i=l 


N 


N - 1 
N 

' N - 1 


(w V Il - ^il v W v I v + 


- T2 w vlu 

W v l v - - 


(52) 

(53) 

(54) 

(55) 


where we use the definition of from Eq. 34 This 


quantity can be calculated directly from the aggregates 
stored in the database. 

Every time a new association is stored in the database, 
the values of the variability parameters calculated at the 
time of association are stored together with it. In this 
way, it is possible to query the database for the variability 
parameters corresponding to any point in the history of a 
particular running catalogue entry. 

6.5. Pixel store 

Images are not created during the operation of the TraP; 
in general, therefore, we regard their storage as outwith the 
scope of TraP operations. However, it is often convenient 
to maintain easy access to image data which has been 
processed. This enables end users who are analysing the 
TraP results to quickly cross-check them with a visual 
inspection of the image data. Indeed, tools such as Banana 


(t 5.2.2) can over-plot details of sources identified by the 


TraP on the image data. 

In normal operation, the TraP reads images from the 
filesystem attached to whichever machine (or machines) 
upon which it is executing. Often, that filesystem is not 
intended as long-term image storage, but is rather a tem¬ 
porary resting place on whatever compute system is being 
used for analysis. Further, it may not always be desir¬ 
able (for security or management reasons) for the ultimate 
scientific user of the TraP to have access to the systems 
upon which the pipeline runs. Finally, it is simply more 
convenient to aggregate images for display in one location, 
rather than have Banana or other tools search for them on 
diverse filesystems. 

For these reasons, the TraP can optionally insert a copy 
of all pixel data it processes to a centralized store. The 
term “pixel data” is used deliberately: rather than stor¬ 
ing complete image cubes, with full metadata, images are 
reduced to a lowest common denominator form consisting 
of just a pixel grid and coordinate system stored in FITS 
format. This enables a convenient and uniform interface 
by which data may be accessed for display, but does not 
amount to a comprehensive archive of the images. 


The pixel storage used by the TraP is implemented as a 
MongoDEp’] database. MongoDB is a “document-oriented” 
database, which makes it easy to simply store and retrieve 
large “blobs” of binary data (such as our pixels) using a 
simple key-value look-up scheme. 

Pixel data may be saved to the MongoDB database by 


the data accessor (14.1I when it is first loaded from disk. A 


URL identifying the location of the corresponding pixels is 
then stored in the Image table of the main TraP database 

(Fig# 

6 .6. Dissemination of transient notifications 

After an event has been selected as scientifically note¬ 
worthy, information about it must be rapidly distributed. 
In general, notifications will be sent to the community at 
large, although it is possible that certain events may only 
be shared with selected partners. 

Currently, the rate of transients being announced by 
LOFAR is low, but we anticipate it increasing in the future 

Look¬ 


(Fender et al. 

in prep.; 

Stewart et al., 

in prep. 


ing further ahead to facilities like SKA (Dewdney et al. 


2010) and LSST (Ivezic et al. 2014), it is reasonable to 


expect that millions of transients may be announced every 
day. Furthermore, rapid turn-around time for follow-up 
observations is often necessary. Therefore, we regard it as 
imperative that, as far as is possible, transient alerts can be 
generated, transmitted, received and acted upon without 
human intervention. This makes possible the development 
of the automatic systems that will be required to handle 


the upcoming transient deluge (see, for example, Staley 


et al. 2013). 


With the above considerations in mind, we have stan¬ 


dardized upon the VOEvent (Seaman et al. 2011) format 


developed by the International Virtual Observatory Al¬ 
liance (IVOA) for describing transients detected by LOFAR. 
VOEvent provides a standardized, machine-readable way of 
describing a celestial event with the implication that timely 
follow-up is of interest. VOEvent provides mechanisms for 
describing: 

• The facility and/or observer responsible for publish¬ 
ing the notification packet; 

• A description of the event observed; 

• Where and when the observations where made; 

• Instrument specific information describing how the 
data was collected; and 

• A scientific assessment of the event, which may be 
used to motivate the request for follow-up. 

All of this information is presented in an XML document 
which can be conveniently manipulated by computer, but 


5 https://www.mongodb.org/ 
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it may also be accompanied by plain text descriptions for 
human consumption. 

The flexibility of this format is such that early LOFAR 
transient notifications can be simple (a position, a times¬ 
tamp, a frequency and a flux density measurement, for 
example), and, as our understanding of both the instru¬ 
mentation and the low frequency radio sky improves, the 
event descriptions can become increasingly elaborate and 
include detailed classification information and scientific 
assessment. 

The VOEvent standard does not specify a means by 
which VOEvents should be transmitted from originator 
to recipient. However, ongoing work in the IVOA and 
the transient astronomy community has developed a trans¬ 
portation protocol (Allan and Denny 2009) and an early 


version of a worldwide distribution network (Williams et al. 
2012). The TKSP team has developed Cornelj^ (Swinbank 


2014), an open-source implementation of this transporta¬ 


tion system, and will use that to publish VOEvents to the 
distribution network. A prototype of a similar system is 
currently being used to notify the robotic pt5m telescopj^] 
of observations by AMI-LA. 


7. Development methods 

The TraP is a large and complex project: it consists of 
some tens of thousands of lines of code, written in Python 
and SQL, which are very different languages; it ingests 
image data from a variety of different sources; it interacts 
with two different types of database; and it is developed, 
tested and supported by a heterogeneous team of soft¬ 
ware developers and academic astronomers spread across 
multiple different institutions. Ensuring the delivery of 
reliable software which produces scientifically valid results 
under these circumstances requires a rigorous development 
methodology. 

7.1. Planning and issue tracking 

TraP releases are made at the cadence of a few per year. 
This provides a compromise between deploying new and 
upgraded features to end users as rapidly as possible, and 
providing a stable base which users can trust to provide 
consistent results from day to day while they work on a 
particular science project. 

Releases alternate between “science” and “technical” 
focuses. The science-focused releases aim to deploy new 
and upgraded scientific analysis capabilities. Technically 
focused releases concentrate on consolidation of the code- 
base and introducing new technology, without changing 
the capabilities available to end users. 


2<: http: //comet .transientskp. org/ 

2 'https://sites.google.com/site/point5metre/ 


Goals for a release are defined through a series of “is¬ 
sues” targeted to a particular milestone in an issue trackeip^] 
During the development cycle a daily test build is made 
available for commissioners. In light of development expe¬ 
rience and results from the test build, the issues targeted 
for the milestone may be revised, and new issues may be 
added. When all the issues targeted for that milestone have 
been addressed, a release occurs and the cycle repeats. 

7.2. Code repository and version control 

The most fundamental tool in developing and maintain¬ 
ing a large codebase is a version control system. This is 
essential to maintain a list of changes to the code, includ¬ 
ing information about who changed what, when, and what 
the rationale was. We use the version control system to 
develop and test multiple variations of the TraP in parallel; 
to isolate and revert errors introduced to the code; and 
to enable the painless integration of code developed by 
different and geographically separate developers. 

The TraP makes use of Gitp 7 ] with a central repository 
currently hosted on Gitl lulf {ll | Our experience has been 
that software developers are quick to adapt to working with 
Git, but that its complexity can be off-putting to those 
coming from a more purely scientific background. We 
have organized training sessions and workshops in order to 
mitigate this. 

7.3. Code review 

In order to ensure that all code entering the codebase is 
of high quality, and to ensure that there is no single part of 
the codebase which is understood by only one developer, we 
require that all code contributions are reviewed by a team 
member other than their original author before they are 
added to the TraP. This process is managed using GitHub’s 
“pull request” interface. 

The overhead introduced by this review step is not 
negligible: the reviewer must often invest considerable 
time to become familiar with the code being reviewed, 
and sometimes a lengthy discussion between the original 
author and the reviewer can result. Furthermore, it can 
occasionally be frustrating for the author to wait for a 
reviewer to become available during busy times. 

Despite these downsides, though, the review process 
has been successful: since it was instituted, the quality and 
reliability of the TraP codebase has increased markedly, 
and the entire development team has better insight into all 
parts of the pipeline rather than just their own particular 
specialization. 


28 For most of the lifetime of the TraP to date, this was Redmine, 
http:///www.redmine.org/ We have recently switched to GitHub 
Issues, http://www.github.com/ for better integration with our ver¬ 
sion control system. 

2! 'http: //www. git-scm. com/ 

* i( ‘http: //www. github. com/ 
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7. ^. Testing and continuous integration 

Testing is fundamental to the development of any soft¬ 
ware system. Mistakes are inevitable, and, in a large and 
complex codebase, predicting all possible effects of even 
simple changes becomes challenging. This is particularly 
the case when development takes place using a dynamic 
language such as Python: with no compile-time checking 
for type or even syntax, it is easy for errors to slip by 
without being noticed. 

The TraP codebase is rigorously and automatically 
tested. At time of writing, the test suite consists of some 
347 individual test cases, with three times that number of 
individual assertions contained within them. Test cases 
cover everything from ‘unit’ testing (checking that individ¬ 
ual functions and procedures perform as expected when 
provided with both normal and extraordinary input condi¬ 
tions) to large scale ‘integration’ tests which validate the 
scientific results produced by large sections of the TraP 
on given input data. All new code must pass all of these 
tests (or, alternatively, explain why the test suite should be 
changed) before it is accepted by the code review process. 
Furthermore, all submissions are expected to come with 
their own set of tests which demonstrate their correctness. 

Our testing infrastructure is based upon the unittest 
module provided as part of Python’s standard library and 
the Jenkin^l continuous integration system. 

7. 5. Documentation 

Documentation is provided both for the end-user as¬ 
tronomer who needs to understand how to process their 
data and interpret the results, and for the expert user 
or developer who is extending the TraP to address their 
particular use case. 

The TraP is documented using the Sphimj^Jdocumenta- 
tion system. This can both automatically generate interface 
documentation from the TraP’s Python code while also in¬ 
corporating hand-written material giving a more complete 
description of the code along with tutorial-style documen¬ 
tation. As part of the code review process, reviewers are 
expected to check that code not only functions properly 
and is well tested, but also that, if appropriate, it is ac¬ 
companied by appropriate additions or alterations to the 
documentation. 

The documentation for all released versions of the TraP 
as well as the latest developmental version is available from 
the project websitep] 

8. Integration testing 

As described in §7.4| the TraP codebase is well covered 
by an extensive test suite which tests individual compo¬ 
nents and their interactions when provided with a variety 


of different inputs, based on both synthetic and archival 
data. Further, individual subsystems and the algorithms 
they implement have undergone extensive testing both in 
the published literature and in regular use. For example, 


Scheers (2011 chapter 4) describes how the source associa¬ 
tion routines were applied to cross-match the VLSS (VLA 


Low Frequency Sky Survey; 

Cohen et al. 

2007), WENSS 

(Westerbork Northern Sky Survey; 

Rengelink et al. 

1997 

and NVSS (NRAO VLA Sky Survey; 

Condon et al. 

1998 


catalogues. This functionality forms the basis of the LO- 


FAR Global Sky Model database ([van Haarlem et al. 2013). 


Similarly, Spreeuw (2010 chapter 3) describes an elaborate 
series of statistical tests on the sourcefinder, which are ex¬ 


panded upon by Carbone et al. (in prep.). Results from the 
sourcefinder were also submitted to the ASKAP /EMFp^j 


Source Finding Data Challenge (Hopkins et al. in prep. 


the final results of this exercise have not yet been published, 
but preliminary indications are that the TraP code has 
performed to a high standard. 

Although the individual components of the TraP are 
well tested, it is useful to consider an integration test, which 
demonstrates the operation of the TraP as a coherent whole 
and provides an indication as to how the results may be 
interpreted. It is stressed that this section serves primarily 
as an illustration of a pipeline run under strictly controlled 
circumstances: we do not attempt to account for complex or 
unexpected behaviour of astronomical sources, as this can 
best be considered by comparing the source behaviour to 
the documented sourcefinder capabilities, database sources 
association behaviour, variability metrics, etc. It is worth 


noting that a companion paper, Rowlinson et al. (in prep. 


expands upon the techniques presented here to establish 
strategies to determine optimal TraP configuration for a 
given dataset given expected source and image character¬ 
istics, while early science results derived from pre-release 


versions of the TraP are now becoming available (Carbone 


et al. submitted). 


8.1. Simulation procedure 

Simulated monochromatic lightcurves representing single¬ 
epoch transients observed at a frequency v were generated. 
Each lightcurve consisted of 20 flux density measurements, 
with the flux density recorded for measurement i, I V} i, given 
by: 

/transient if i — 8, / -\ 

I quiescent otherwise. 

The transient flux density, /transient was varied over the 
range [5, 95] Jy in steps of 5 Jy. The quiescent flux density, 
/quiescent, was varied over the range [0 Jy, /transient) using 
the same step size. In this way, a total of 190 lightcurves 
were generated. 

For each lightcurve, a set of 20 images representing 
LOFAR observations of the transient was simulated. In 


31 http://jenkins-ci.org/ _ 

32 http://sphinx-doc.org/ 34 EMU is ASKAP’s Evolutionary Map of the Universe Survey 

3J http://docs.transientskp.org/ Science Project 
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order to closely mimic genuine LOFAR observations, the 
simulation developed was based closely on the structure of 
existing LOFAR visibility data. We started with visibility 
data obtained as part of LOFAR’s Multifrequency Snapshot 


Sky Survey (MSSS; Heald et al. 20141. The data consisted 
of 20 observations of the field of 3C295 (14 h, ll m 20.6 s , 
+52°12 , 21") made between 03:00Z and 08:00Z on December 
24, 2011. Each observation had an integration time of 11 
minutes, and was followed by a 4 minute re-pointing time. 
Observation configurations were identical, covering 2 MHz 
of bandwidth, divided among 10 subbands and centred on 
54 MHz. 

We generated a model sky for the area being observed by 


selecting all sources from the VLSS catalogue (Cohen et al. 


2007) which fall within 8 ° of the pointing centre and which 


are above a limiting flux density of 1 Jy. Spectral indices 
for these sources were generated by comparing the VLSS 


flux densities with those reported in WENSS (Rengelink 


et al. 1997) and NYSS (Condon et al. 1998), then used 


to extrapolate the source flux density to 54 MHz. Where a 
VLSS source had no counterpart in the other catalogues, a 
spectral index of -0.7 was assumed. 

The simulation procedure for each image was: 

All subbands were averaged to produce a single chan¬ 
nel with a width of 1 MHz; 

The stored visibilities were replaced with randomly 
generated Gaussian noise at a level chosen to match 
the System Equivalent Flux Density of the instrument 


1 . 


2 . 


(van Haarlem et al. 2013); 


3. 


The appropriate transient flux density, based on the 
lightcurve being processed and the image number, was 
appended to the model sky at position 14 /l 20 m 00.0 s , 
+52°00'00.0" (1.34° from 3C295); 

4. BBS (Loose} [2008|, the standard tool use for cali- 


5. 


brating LOFAR data, was used to simulate model 
visibilities and add them to the data based on a 
user-supplied model sky; 

The data was calibrated and imaged as usual. The 
transient source was included in the model sky used 
for calibration. Images were generated with a radius 
of 6 °. 


When calculating 77 ^ (Eq. 36), we assigned an equal 
weight (equal to the reciprocal of the average error 
across all flux density measurements) to each data 
point. 


• Only those flux density measurements recorded in 
or after the image of first detection are included in 
the variability metric calculation; as per (|3j measure¬ 
ments from earlier images are not available during 
pipeline processing. 


In this way, we were able to predict the variability 
metrics which the TraP should calculate for each source, 
and determine in advance which ones ought to be identified 
as transients for a given TraP configuration. 


8.3. Transients pipeline results 

Each dataset of snapshot images (190 datasets, 1 for 
each transient) was run through TraP. A near-default 
pipeline configuration was used: the quality control system 
set to not reject any images and the shape of all point 
sources was constrained to be equal to that of the clean 
beam. The variability metrics corresponding to the final 
snapshot for each source were extracted from the database 
at the end of each pipeline run. Fig. [ 6 ] shows the predicted 
values for each of the transient sources in comparison to 
the measured value by the transient pipeline. The scatter 
in this figure is due to three factors: 


• The simulation process generates sources on a noisy 
background, and this noise impacts on the results 
produced by the sourcefinder; 


When a source is not detected by the initial sourcefind¬ 
ing step, the TraP’s null-detection procedure ((4.5) 
will force a constrained fit to its position and record 
a (likely non-zero) flux density. The prediction pro¬ 
cedure, by contrast, assumes a flux density of 0 Jy; 


• The predicted values were calculated assuming an 
equal weighting for all flux density points, whereas 
the TraP assigns each point an independent weight. 


The resulting images had an RMS noise level around 
0.5Jy/beam. Simulated sources with /quiescent > 5 are 
therefore detected at around 10 cr or higher in their quiescent 
states; sources with / qu i esce nt = 0 are detected when the 
transient turns on in image 8 . 

8.2. Predicted results 

The flux density coefficient of variation, V v , and the 
significance of the variability, 77,,, as described in §4.7.2| were 
calculated independently of the pipeline machinery for each 
of the 190 lightcurves described in the previous section. 
Note that the calculation of these metrics depends not only 
on the raw simulated values, as described in the previous 
section, but also on the configuration of the pipeline run. 
In particular: 


Taking these expected deviations into account, the pre¬ 
dicted and measured values are consistent with each other. 
In particular, we tend to predict higher variability indices 
than are measured for faint sources, as our prediction pro¬ 
cedure takes no account of image noise, assigns an equal 
weight to all measurements, and assigns a flux density of 
0 Jy to non-detections. 

As per (5.1 we distinguish between two classes of tran¬ 
sients: new detections and variable sources. In this section 
we describe targeted tests which confirm that both mecha¬ 
nisms are performing as expected. 

Throughout this section, we refer to “true positive” (TP) 
detections when a transient source is correctly identified as 
such; “false positive” (FP) detections when a non-transient 
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(a) Flux density coefficient of variation (14). 



Figure 6: Comparison of predicted and measured variability parameters for simulated transients. 


source is incorrectly labelled as transient, and “false nega¬ 
tive” (FN) detections when a transient source is incorrectly 
labelled as non-transient. We then define 


Precision = 
Recall = 


TP 

TP + FP 
TP 

TP + FN' 


and 


(57) 

(58) 


Following these definitions, the precision of the result 
is the fraction of the total number of detections which are 
correct and the recall is the fraction of the total number 
of transients which were correctly identified. The best 
possible TraP performance is obtained by maximizing both 
the precision and the recall. 


8.3.1. New detections 

We selected all sources from the database that were 
initially detected in any image except the first. Following 
the procedure described in (5.1.1 and using a margin of 
3<r, we classified them as either not transient (i.e. below 
previous detection threshold), marginally transient (above 
previous detection threshold in the lowest-noise portion 
of at least one previous image) or likely transient (above 
previous detection threshold in the highest-noise portion 
of at least one previous image). We find that: 


• Likely transients are recovered with a precision of 1.0 
and recall of 0.94; 

• Marginal transients are recovered with a precision of 
0.02 and recall of 1.0. 


Based on the simplified test described here, we conclude 
that the algorithm used to detect likely transients provides 
a robust way of identifying many transients with a high re¬ 
sistance to false positives. Further, since all lightcurve data 
is retained in the database, the list of marginal transients 
provides a key starting point for future manual checking 
and data mining. 


An important limiting factor in this test is the limited 
resolution at which the noise maps are stored (i.e. just a 
“best” and “worst” value for each image). Recording noise 
at a more fine-grained level would enable us to significantly 
increase the precision with which possible transients are 
identified. This is a possible area of future TraP develop¬ 
ment (p|. 

8.3.2. Variable sources 

We selected all sources which were initially observed 
in the first image (i.e. they were not candidates for being 
marked as new detections) and which had values of rp, and 
V v greater than 0.1 and 0.01 respectively. We constructed 
an equivalent list based on the simulation inputs and known 
image noise level; note that this list excludes some transient 
sources which fall below the detection threshold. By com¬ 
bining these lists, we can calculate the precision and recall 


(Eqs. 57 & 58) as a function of the variability parameters. 
These are plotted in Fig. [7} Above some limiting value of 
each threshold, there are no positive detections (either true 


or false) so the values of Eqs. and 58 are undefined; the 


plots are truncated at this point. Note that for values of 
? 7 „ > 1 and V v > 0.3 the precision is 1.0: no false positives 
are recorded. Below these values, precision drops rapidly 
due to noise-based variation of stable sources. 

For all values of r\ v and V„ the recall is similarly close to 
1.0. Variations are due to uncertainties introduced by the 
simulation and measurement process, which occasionally 
cause the measured value of the transient parameters to 
drop below their predicted values. 


9. Run-time performance 

As described in sections [272] and [3| the TraP is ultimately 
intended to perform near real-time analysis of streaming 
image data. Although the required rapid imaging capabil¬ 
ity is not yet available from LOFAR, we anticipate that 


other projects—most notably AARTFAAC (Prasad and 


Wijnholds 2012) -will provide streaming image data in 
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(a) Probabilities as a function of the threshold applied to the 
flux density coefficient of variation (V"„). 


Figure 7: The precision and recall probabilities as functions 
nv- 



(b) Probabilities as a function of the threshold applied to the 
significance of the flux density variability (r/ v ). 

the triggering thresholds for the variability metrics V v and 


the relatively near future. It is therefore instructive to 
consider to what extent the run-time performance of the 
TraP is adequate to address such a data stream. 

It is worth emphasizing that TraP development to date 
has focused on correctness rather than performance. Our 
aim has been to produce a robust and well-tested code- 
base that can then be optimized to address real-time data 
processing. We strive to adopt fundamental algorithms 
which show benign scaling characteristics to large numbers 
of images and sources, but emphasize that the codebase 
still provides ample opportunity for optimization. 

Broadly, we split our consideration of performance char¬ 
acteristics into two parts, corresponding to the two most 
computational expensive parts of the TraP. In §9.1| we 
focus on the performance of the Python code, and particu¬ 
larly that of the computationally intensive source finding 
algorithms. In §9.2| we turn to the matter of inserting 
and associating measurements in the database and the 
calculation of per-source aggregates. 

Throughout, we emphasize that there are a large num¬ 
ber of tunable parameters in this analysis, both in terms 
of the pipeline configuration and the characteristics of the 
test data; here, we only give an overview of likely scaling 
considerations. For a detailed review of the sourcefindcr 
performance refer to Carbone et al. (in prep.), and for an 
in-depth discussion of database characteristics see |Scheers| 


et al. (in prep. 


Throughout this section, the times reported correspond 
to Python code running on an Intel Xeon E5-2660v2 CPU 
with a maximum clock speed of 2.2 GHz and, where ap¬ 
plicable, interoperating with PostgreSQL 9.3.5 running on 
an AMD Opteron 2384 with a maximum clock speed of 
2.7 GHz. We configured PostgreSQL to make better use 
of the available system resources by increasing its working 
memory (to 100MB), its shared buffer (to 2048MB) and 
its checkpoint interval to 32 segments. 



Number of pixels 


Figure 8: Time taken to calculate the background and 
RMS maps as a function of the number of pixels in the 
image. The solid line shows the mean measured time over 
3800 test images; the dashed line, a quadratic fit to the 
data. 


9.1. Sourcefinder performance 

Based on the discussion in §4.3[ we divide the opera¬ 
tion of the sourcefinder into two major components: the 
calculation of per-image background and RMS maps, then 
the identification and measurement of sources within the 
image. The former depends on the size of the image, but is 
independent of the number of sources within it; the latter 
increases with source count. 

For each of the 3800 images simulated as described 
in §8.1| the time taken to generate background and RMS 
maps covering the whole image was measured. The edges 
of the images were then masked, and timing repeated for 
maps covering only the unmasked portion. This process 
was repeated until only a small fraction of the image was 
left unmasked. The times recorded are shown in Fig. [8j 
For comparison, we also plot the results of a least squares 
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Figure 9: Time taken to find and measure the sources in an 
image as a function of the number of sources found. The 
solid line shows the mean time measured over 3800 test 
images; the dashed line a linear fit to the data. 


quadratic fit to the data: 

f map = 1.9 x 10" 1 V - 2.5 x 10 ~ 6 p + 0.1 s, (59) 

where f map is the time taken to process p pixels. While 
the detailed values are system dependent, it is important 
to note that the algorithm scales as 0(N 2 ) in number of 
pixels. 

As the unmasked area of each image is decreased, the 
number of sources which can be detected and measured 
within the image also decreases. For each portion of each 
image we performed a source finding and measurement step 
with a detection threshold of 10er. The sourcefinder was 


configured not to deblend sources (14.3.31, and to constrain 
the shape of the resulting measurements to be equal to the 
restoring beam (j |4.3.5 ). 

Fig. [9] records the total time taken to identify and fit all 
the sources in an image as a function of the source number. 
A linear least squares fit to the data provides the expression 


ta t = 0.012n +0.053 s 


(60) 


for the ifit taken to identify and measure n sources. While 
again the detailed timings are system-dependent, the key 
point is the scaling as O(N) with number of sources. 


9.2. Database performance 

There are two important axes along which database 
performance could vary. The first is with number of images 
processed: as more data is stored, the number of source 
measurements which must be associated and the number 
of data points of which aggregates must be calculated 
increases. For use in a long term monitoring programme, 
we require that this accumulation of data does not cause the 
database to become slower with time. Secondly, we consider 
performance as a function of the number of sources per 
image: more measurements increase not only the number of 


Figure 10: Time taken to process each simulated observa¬ 
tion with a population of 1100 sources in the database. 


aggregates to be calculated but also the number of potential 
source associations. 

Artificial source lists representing an artificial sky at 
arbitrary frequency and pointing and covering a circular 
region of radius 20° were constructed. Sources were placed 
on a regular grid within the region. Each source was 
assumed to be a point source, and assigned a random flux 
density in the range 1-10 Jy from a uniform distribution. 
Sixteen such source lists were created, containing between 
50 and 1200 sources in increments of 50. 

For each source population, a set of 100 source measure¬ 
ments was constructed by perturbing the source position 


with a Fisher distribution (Fisher. 1953) with concentra¬ 
tion parameter k = 2x 10 , chosen to approximate the 
systematic position uncertainty of around 5” which we 


have observed to be typical in LOFAR images (Carbone 


et al. 


m prep. 


. This simulates an observation of the source 

population. 

For each source population, each simulated observa¬ 
tion in turn was inserted and processed (including source 
association and calculation of aggregate parameters by 
the database). The time taken to perform all database 
operations was recorded. 

In Fig. [lO] we show the time taken to perform all the 
processing of each simulated observation of 1100 sources as 
a function of image number; a similar pattern is observed 
for all other source counts. The characteristic “saw tooth” 
pattern in the figure is due to PostgreSQL periodically 
checkpointing its write ahead log; other minor variations 
are explained by internal housekeeping tasks running within 
the database and by varying system and network load over 
the course of the test. The key result, though, is that there 
is no evidence of a systematic increase in processing time 
with observation number. 

In Fig. we show the mean time taken to perform all 
the processing of each simulated observation as a function 
of source count. For comparison, we also plot the expression 


t dh = 0.0017n + 0.10 s 


(61) 
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Number of sources 

Figure 11: Time taken by the database to process a sim¬ 
ulated observation as a function of number of sources in 
the observation. The solid line shows the mean value over 
100 test observations; the dashed line, a linear function for 
comparison. 


where fdb is the time taken to process a simulated observa¬ 
tion of n sources. The detailed timings are, again, system 
dependent, but it is important to note the scaling as O(N) 
with number of sources per observation. 


9.3. Practical performance considerations 

We conclude this discussion of pipeline runtime perfor¬ 
mance by comparing the measured TraP performance to 
potential LOFAR transient monitoring strategies. 


The initial Radio Sky Monitor strategy (Fender et al. 


m prep. 


is to use six beams from LOFAR to tile out a 
wide area on the sky. Each beam consists of four frequency 
bands, which are imaged separately. Assuming LOFAR 


is operating in 8 bit mode (12.1), each band contains 20 


subbands and provides a bandwidth of 3.6 MHz. 

In Tab. [2] we list the parameters of each of six major 
LOFAR observing modes and provide the full width at 
half maximum and the angular resolution of a single image 
constructed using the survey strategy described. The ob¬ 
serving modes include using only the core LOFAR stations, 
using the full Dutch LOFAR array, and using only that 
subset of the full array which contains baselines no more 
than 6 km in length. This latter configuration has been 
shown to provide a good compromise between image fidelity 
and processing time in early LOFAR observations and has 
been used for initial RSM observations. 

In Tab.[3]we provide estimates of the 5er detection limit 
and corresponding source count for each of the configu¬ 
rations at a range of integration times. Sensitivity and 
confusion limits were estimated using the online LOFAR 
Image Noise Calculatoip’l Following current standard LO¬ 
FAR observing practice, when using the high band remote 
stations were “tapered”, reducing sensitivity while increas¬ 
ing field of view to match that of the stations in the core. 


35 Version 0.31; http://www.astron.nl/~heald/test/sens2.php 


Where possible, source counts were taken from preliminary 


analysis of LOFAR RSM data (Fender et al. in prep. 


otherwise, they were estimated by extrapolating from other 
surveys assuming a spectral index of —0.7: for each configu¬ 
ration, the survey best approximating the spatial resolution 
was selected from FIRST (Faint Images of the Radio Sky 
at Twenty-cm; Becker et al. 1994| ), NVSS (Condon et al. 


l998j) , VLSS ( jCohen et al.H2007| ) and WENSS ( jRengelink 


et al. 1997). 


Also in Tab. [3] we estimate the processing time for each 
image through both the sourcefinder (t map + ifit) and the 
database (fdb)- Note that these equations reflect the par¬ 
ticular characteristics of the systems used for benchmark¬ 
ing and can provide only guideline performance estimates. 
Further note that this estimate includes only the major 
pipeline components identified above; it is reasonable to 
expect modest overheads from other parts of the pipeline. 

Note that for most image types excluding the full array 
the processing time per image is considerably less than 
the integration time for that image. This indicates that 
processing a real-time stream of images of this type is 
tractable on the systems used for benchmarking. 

As described above, the likely operation configuration 
consists not of a single image stream but of four bands in 
each of six beams: a factor of 24 greater than the figures 
quoted. For all 1 s, and some longer, integrations, this 
increases the total processing time to greater than the 
integration time. However, three factors mitigate this: 


• As described in j ]6.3| we support parallel (and, option¬ 
ally, distributed) sourcefinder operation: by distribut¬ 
ing the processing time required over multiple CPUs, 
we achieve a near-linear reduction in wall-clock time. 

• By starting to search for sources in timestep t n+ i 
before timestep t n has finished processing, we can 
ensure that both sourcefinding and database systems 
are efficiently occupied, rather than synchronously 
waiting for each other. 

• The hardware platform used for benchmarking the 
database is several years old; more powerful systems 
are now available which will provide a significant 
constant factor improvement in database throughput. 


We also emphasize that the real-time LOFAR imaging 
mode is still under development ((2.2): all current process¬ 
ing is carried out off-line. The run-time performance of the 
TraP is significantly better than that of the imaging tools 
used to provide it with input in this mode. 

Finally, we reiterate that TraP development to date has 
focused on correctness over performance. There is much 
still to be gained in terms of optimization of individual 
algorithms, efficiency of implementation (e.g. replacing core 
Python loops with Cython equivalents) and best exploiting 
the performance of high level tools (e.g. the potential gains 
of the MonetDB system over PostgreSQL as described in 

pXT). 


26 

































Configuration 

name 

Freq. 

(MHz) 

BW 

(MHz) 

Num. 

Core 

stations 

Remote 

Max. baseline 
(km) 

FWHM 

(deg.) 

Ang. res. 
(asec.) 

Pixels per 
image 

LBA core only 

60 

3.6 

24 

0 

3.5 

9.8 

294.7 

1.0 x 10 5 

LB A 6 km baselines 

60 

3.6 

24 

4 

6.0 

9.8 

171.9 

3.0 x 10 5 

LBA NL array 

60 

3.6 

24 

16 

121.0 

9.8 

8.5 

1.2 x 10 s 

HBA core only 

150 

3.6 

24 

0 

3.5 

3.8 

117.9 

9.5 x 10 4 

HBA 6 km baselines 

150 

3.6 

24 

4 

6.0 

3.8 

68.8 

2.8 x 10 5 

HBA NL array 

150 

3.6 

24 

16 

121.0 

3.8 

3.4 

1.1 x 10 8 


Table 2: Parameters of LOFAR observing modes. For each mode, we quote the angular resolution and full width at half 
maximum of a single beam and a bandwidth of 3.6 MHz, equivalent to 20 subbands. 


Configuration 

Integration 

5cr detection limit 

Source count 

Processing time (s) 

name 

time (s) 

(mjy/beam) 

Origin 

Number 

Sourcehnder 

Database 

LBA core only 

1 

2346.5 

VLSS 

27.6 

0.46 

0.15 


10 

C 1066.3 

VLSS 

90.0 

1.22 

0.25 

LBA 6 km baselines 

1 

2005.1 

VLSS 

34.9 

1.58 

0.16 


10 

634.1 

VLSS 

196.3 

3.55 

0.42 


100 

C 465.0 

VLSS 

312.6 

4.98 

0.62 

LBA NL array 

1 

1395.8 

VLSS 

60.1 

2.86 x 10 5 

0.20 


10 

441.4 

VLSS 

337.9 

2.86 x 10 5 

0.67 


100 

139.6 

RSM 

761.3 

2.86 x 10 5 

1.37 


1000 

44.2 

WENSS 

2043.6 

2.86 x 10 5 

3.50 


10000 

14.0 

FIRST 

3537.0 

2.86 x 10 5 

5.99 

HBA core only 

1 

C 137.0 

LOFAR 

45.3 

0.67 

0.18 

HBA 6 km baselines 

1 

102.0 

LOFAR 

70.5 

1.86 

0.22 


10 

C 59.7 

LOFAR 

157.4 

2.92 

0.37 

HBA NL array 

1 

82.8 

LOFAR 

96.4 

2.53 x 10 5 

0.26 


10 

26.2 

WENSS 

323.0 

2.53 x 10 5 

0.64 


100 

8.3 

FIRST 

449.5 

2.53 x 10 5 

0.85 


1000 

2.6 

FIRST 

2516.8 

2.53 x 10 5 

4.29 


10000 

0.8 

FIRST 

14153.5 

2.53 x 10 5 

23.64 


Table 3: Source counts and processing times predicted for a single image using each of the LOFAR configurations 
described in Tab. [2] at a range of integration times. The symbol C indicates that the detection limit was set by confusion 
noise rather than image sensitivity. 
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10. Releases and future development 

This manuscript describes release 2.0 of the TraP, dating 
from December 2014. With this release, the TraP became 
an open source project, and therefore freely available to 
download from our GitHub repositorjj^] under a two-clause 
BSD-stykFH license. If you use the TraP, or code derived 
from it, in a paper or other publication, we request that 
you cite this work. 

After this release, development will continue. Future 
releases of the TraP are expected to expand upon the 
current functionality to offer features such as: 

• Performance optimization in support of real-time 
streaming transient monitoring; 

• Automatic preliminary classification of detected tran¬ 
sients; 

• Support for multi-terabyte lightcurve archives and 
enhanced archive data-mining and visualization; 

• Automatic cross-correlation of TraP detected sources 
with known catalogue sources across a range of wave¬ 
lengths; 


transients for offline analysis. We have demonstrated that 
it is capable of accurately recovering a known population 
of transients from simulated LOFAR observations. 

The TraP is now being used in support of ongoing 
LOFAR observing campaigns. However, development con¬ 
tinues, and we are actively expanding its capabilities, both 
to better address ongoing LOFAR operations and to in¬ 
crease its applicability to other instruments and wavelength 
regimes. The codebase is open source and freely available; 
we actively invite you to join us in improving it. 
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• Improved methods for apportioning flux density from 
a single measurement among multiple sources on 
association; 

• Full Stokes support, for both identifying and classify¬ 
ing sources; 

• Higher resolution and more flexible noise maps; 

• Construction and usage of a deep average image of 
the surveyed area. 

11. Conclusions 

The current and next generation of astronomical survey 
facilities, across a wide range of wavelengths but particu¬ 
larly in the radio, provide an opportunity to explore the 
transient and variable sky in powerful and unprecedented 
ways. This is especially true of LOFAR, which combines a 
remarkably wide field of view with unique sensitivity to low 
radio frequencies and a flexible, software-driven architec¬ 
ture. However, identifying transients in the massive data 
volumes produced by these instruments is challenging. 

This manuscript has described our attempt to rise to 
this challenge in the form of the LOFAR Transients Pipeline. 
It combines a flexible, high-performance architecture with 
robust analysis tools in a well tested and documented 
package. We have shown how it can both be used to 
generate alert messages as new transients are discovered 
and to populate a database of lightcurves of potential 
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